Changeset 7276


Ignore:
Timestamp:
Jun 30, 2009, 2:07:41 PM (14 years ago)
Author:
ole
Message:

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

Files:
1 deleted
201 edited
181 copied

Legend:

Unmodified
Added
Removed

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

    r7053 r7276  
    2626from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    2727     import Transmissive_boundary
    28 
    2928from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain
    3029from anuga.abstract_2d_finite_volumes.region\
    … …  
    3433from quantity import Quantity
    3534
    36 import Numeric as num
     35import numpy as num
    3736
    3837
    … …  
    9594        """
    9695
    97 
    9896        number_of_full_nodes=None
    9997        number_of_full_triangles=None
    … …  
    193191            buffer_shape = self.full_send_dict[key][0].shape[0]
    194192            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
    195                                                       num.Float))
     193                                                      num.float))
    196194
    197195        for key in self.ghost_recv_dict:
    198196            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
    199197            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys),
    200                                              num.Float))
     198                                             num.float))
    201199
    202200        # Setup cell full flag
    … …  
    205203        N = len(self) #number_of_elements
    206204        self.number_of_elements = N
    207         self.tri_full_flag = num.ones(N, num.Int)
     205        self.tri_full_flag = num.ones(N, num.int)
    208206        for i in self.ghost_recv_dict.keys():
    209207            for id in self.ghost_recv_dict[i][0]:
    … …  
    267265        # (boolean) array, to be used during the flux calculation.
    268266        N = len(self) # Number_of_triangles
    269         self.already_computed_flux = num.zeros((N, 3), num.Int)
     267        self.already_computed_flux = num.zeros((N, 3), num.int)
    270268
    271269        # Storage for maximal speeds computed for each triangle by
    272270        # compute_fluxes.
    273271        # This is used for diagnostics only (reset at every yieldstep)
    274         self.max_speed = num.zeros(N, num.Float)
     272        self.max_speed = num.zeros(N, num.float)
    275273
    276274        if mesh_filename is not None:
    … …  
    391389            raise Exception, msg
    392390
    393         q = num.zeros(len(self.conserved_quantities), num.Float)
     391        q = num.zeros(len(self.conserved_quantities), num.float)
    394392
    395393        for i, name in enumerate(self.conserved_quantities):
    … …  
    463461
    464462        name:  Name of quantity
    465         value: Compatible list, Numeric array, const or function (see below)
     463        value: Compatible list, numeric array, const or function (see below)
    466464
    467465        The values will be stored in elements following their internal ordering.
    … …  
    895893
    896894            msg += '------------------------------------------------\n'
    897             msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed),
    898                                                max(self.max_speed))
     895            msg += '  Speeds in [%f, %f]\n' % (num.min(self.max_speed),
     896                                               num.max(self.max_speed))
    899897            msg += '  Histogram:\n'
    900898
    … …  
    908906                else:
    909907                    # Closed upper interval
    910                     hi = max(self.max_speed)
     908                    hi = num.max(self.max_speed)
    911909                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
    912910
    913             N = len(self.max_speed)
     911            N = len(self.max_speed.flat)
    914912            if N > 10:
    915913                msg += '  Percentiles (10%):\n'
    … …  
    13731371                self.number_of_steps = 0
    13741372                self.number_of_first_order_steps = 0
    1375                 self.max_speed = num.zeros(N, num.Float)
     1373                self.max_speed = num.zeros(N, num.float)
    13761374
    13771375    ##
    … …  
    16421640        # momentum
    16431641        if self.protect_against_isolated_degenerate_timesteps is True and \
    1644                self.max_speed > 10.0: # FIXME (Ole): Make this configurable
     1642               num.max(self.max_speed) > 10.0: # FIXME (Ole): Make this configurable
    16451643
    16461644            # Setup 10 bins for speed histogram
    … …  
    16571655
    16581656                    # Find triangles in last bin
    1659                     # FIXME - speed up using Numeric
     1657                    # FIXME - speed up using numeric package
    16601658                    d = 0
    16611659                    for i in range(self.number_of_full_triangles):
    … …  
    17631761        # We have a list with ghosts expecting updates
    17641762
    1765 
    1766         from Numeric import take,put
    1767 
    1768 
    17691763        #Update of ghost cells
    17701764        iproc = self.processor
    … …  
    17791773            for i, q in enumerate(self.conserved_quantities):
    17801774                Q_cv =  self.quantities[q].centroid_values
    1781                 put(Q_cv,     Idg, take(Q_cv,     Idf))
     1775                num.put(Q_cv, Idg, num.take(Q_cv, Idf, axis=0))
    17821776
    17831777 

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

    r6161 r7276  
    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

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

    r6654 r7276  
    1 import Numeric as num
     1import copy
     2import numpy as num
    23
    34from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    1011    which is represented as a coordinate set (x,y).
    1112    Vertices from different triangles can point to the same node.
    12     The nodes are implemented as an Nx2 Numeric array containing the
     13    The nodes are implemented as an Nx2 numeric array containing the
    1314    x and y coordinates.
    1415
    … …  
    1920    where
    2021
    21       nodes is either a list of 2-tuples or an Nx2 Numeric array of
     22      nodes is either a list of 2-tuples or an Nx2 numeric array of
    2223      floats representing all x, y coordinates in the mesh.
    2324
    24       triangles is either a list of 3-tuples or an Mx3 Numeric array of
     25      triangles is either a list of 3-tuples or an Mx3 numeric array of
    2526      integers representing indices of all vertices in the mesh.
    2627      Each vertex is identified by its index i in [0, N-1].
    … …  
    3738        triangles = [ [1,0,2], [1,2,3] ]   # bac, bce
    3839
    39         # Create mesh with two triangles: bac and bce   
     40        # Create mesh with two triangles: bac and bce
    4041        mesh = Mesh(nodes, triangles)
    4142
    … …  
    5657    """
    5758
    58     #FIXME: It would be a good idea to use geospatial data as an alternative
    59     #input
    60     def __init__(self, nodes, triangles,
    61                  geo_reference=None,                 
     59    # FIXME: It would be a good idea to use geospatial data as an alternative
     60    #        input
     61    def __init__(self,
     62                 nodes,
     63                 triangles,
     64                 geo_reference=None,
    6265                 number_of_full_nodes=None,
    63                  number_of_full_triangles=None,                 
     66                 number_of_full_triangles=None,
    6467                 verbose=False):
    6568        """Build triangular 2d mesh from nodes and triangle information
    6669
    6770        Input:
    68        
     71
    6972          nodes: x,y coordinates represented as a sequence of 2-tuples or
    70                  a Nx2 Numeric array of floats.
    71                  
    72           triangles: sequence of 3-tuples or Mx3 Numeric array of
     73                 a Nx2 numeric array of floats.
     74
     75          triangles: sequence of 3-tuples or Mx3 numeric array of
    7376                     non-negative integers representing indices into
    7477                     the nodes array.
    75        
     78
    7679          georeference (optional): If specified coordinates are
    7780          assumed to be relative to this origin.
    … …  
    8386        In this case it is usefull to specify the number of real (called full)
    8487        nodes and triangles. If omitted they will default to all.
    85          
    86         """
    87 
    88         if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
    89 
    90         self.triangles = num.array(triangles, num.Int)
    91         self.nodes = num.array(nodes, num.Float)
    92 
    93 
    94         # Register number of elements and nodes
     88
     89        """
     90
     91        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
     92
     93        self.triangles = num.array(triangles, num.int)
     94        self.nodes = num.array(nodes, num.float)
     95
     96        # Register number of elements and nodes
    9597        self.number_of_triangles = N = self.triangles.shape[0]
    96         self.number_of_nodes = self.nodes.shape[0]       
    97 
    98        
     98        self.number_of_nodes = self.nodes.shape[0]
    9999
    100100        if number_of_full_nodes is None:
    … …  
    102102        else:
    103103            assert int(number_of_full_nodes)
    104             self.number_of_full_nodes = number_of_full_nodes           
    105 
     104            self.number_of_full_nodes = number_of_full_nodes
    106105
    107106        if number_of_full_triangles is None:
    108             self.number_of_full_triangles = self.number_of_triangles           
     107            self.number_of_full_triangles = self.number_of_triangles
    109108        else:
    110             assert int(number_of_full_triangles)           
     109            assert int(number_of_full_triangles)
    111110            self.number_of_full_triangles = number_of_full_triangles
    112        
    113        
    114         #print self.number_of_full_nodes, self.number_of_nodes
    115         #print self.number_of_full_triangles, self.number_of_triangles
    116        
    117            
    118111
    119112        # FIXME: this stores a geo_reference, but when coords are returned
    120113        # This geo_ref is not taken into account!
    121114        if geo_reference is None:
    122             self.geo_reference = Geo_reference() #Use defaults
     115            self.geo_reference = Geo_reference()    # Use defaults
    123116        else:
    124117            self.geo_reference = geo_reference
    125118
    126119        # Input checks
    127         msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
    128         msg += 'The supplied array has the shape: %s'\
    129                %str(self.triangles.shape)
     120        msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
     121               'The supplied array has the shape: %s'
     122               % str(self.triangles.shape))
    130123        assert len(self.triangles.shape) == 2, msg
    131124
    132         msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
    133         msg += 'The supplied array has the shape: %s'\
    134                %str(self.nodes.shape)
     125        msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples'
     126               'The supplied array has the shape: %s' % str(self.nodes.shape))
    135127        assert len(self.nodes.shape) == 2, msg
    136128
    137129        msg = 'Vertex indices reference non-existing coordinate sets'
    138         assert max(self.triangles.flat) < self.nodes.shape[0], msg
    139 
     130        assert num.max(self.triangles) < self.nodes.shape[0], msg
    140131
    141132        # FIXME: Maybe move to statistics?
    142133        # Or use with get_extent
    143         xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]) ,
    144                       max(self.nodes[:,0]), max(self.nodes[:,1]) ]
    145 
    146         self.xy_extent = num.array(xy_extent, num.Float)
    147 
     134        xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]),
     135                     max(self.nodes[:,0]), max(self.nodes[:,1])]
     136
     137        self.xy_extent = num.array(xy_extent, num.float)
    148138
    149139        # 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)
     140        self.normals = num.zeros((N, 6), num.float)
     141        self.areas = num.zeros(N, num.float)
     142        self.edgelengths = num.zeros((N, 3), num.float)
    153143
    154144        # Get x,y coordinates for all triangles and store
    155145        self.vertex_coordinates = V = self.compute_vertex_coordinates()
    156146
    157 
    158147        # Initialise each triangle
    159148        if verbose:
    160             print 'General_mesh: Computing areas, normals and edgelenghts'
    161            
     149            print 'General_mesh: Computing areas, normals and edgelengths'
     150
    162151        for i in range(N):
    163             if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
    164            
     152            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N)
     153
    165154            x0, y0 = V[3*i, :]
    166155            x1, y1 = V[3*i+1, :]
    167             x2, y2 = V[3*i+2, :]           
     156            x2, y2 = V[3*i+2, :]
    168157
    169158            # Area
    170             self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
    171 
    172             msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
    173             msg += ' is degenerate:  area == %f' %self.areas[i]
     159            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
     160
     161            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2)
     162            msg += ' is degenerate:  area == %f' % self.areas[i]
    174163            assert self.areas[i] > 0.0, msg
    175 
    176164
    177165            # Normals
    … …  
    184172            #     the first vertex, etc)
    185173            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
    186 
    187             n0 = num.array([x2 - x1, y2 - y1], num.Float)
     174            n0 = num.array([x2-x1, y2-y1], num.float)
    188175            l0 = num.sqrt(num.sum(n0**2))
    189176
    190             n1 = num.array([x0 - x2, y0 - y2], num.Float)
     177            n1 = num.array([x0-x2, y0-y2], num.float)
    191178            l1 = num.sqrt(num.sum(n1**2))
    192179
    193             n2 = num.array([x1 - x0, y1 - y0], num.Float)
     180            n2 = num.array([x1-x0, y1-y0], num.float)
    194181            l2 = num.sqrt(num.sum(n2**2))
    195182
    … …  
    207194            self.edgelengths[i, :] = [l0, l1, l2]
    208195
    209        
    210         # Build structure listing which trianglse belong to which node.
    211         if verbose: print 'Building inverted triangle structure'         
     196        # Build structure listing which triangles belong to which node.
     197        if verbose: print 'Building inverted triangle structure'
    212198        self.build_inverted_triangle_structure()
    213 
    214            
    215199
    216200    def __len__(self):
    217201        return self.number_of_triangles
    218    
    219202
    220203    def __repr__(self):
    221         return 'Mesh: %d vertices, %d triangles'\
    222                %(self.nodes.shape[0], len(self))
     204        return ('Mesh: %d vertices, %d triangles'
     205                % (self.nodes.shape[0], len(self)))
    223206
    224207    def get_normals(self):
    225208        """Return all normal vectors.
     209
    226210        Return normal vectors for all triangles as an Nx6 array
    227211        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
    228212        """
     213
    229214        return self.normals
    230 
    231215
    232216    def get_normal(self, i, j):
    233217        """Return normal vector j of the i'th triangle.
     218
    234219        Return value is the numeric array slice [x, y]
    235220        """
     221
    236222        return self.normals[i, 2*j:2*j+2]
    237223       
    … …  
    245231    def get_number_of_nodes(self):
    246232        return self.number_of_nodes
    247        
     233
    248234    def get_nodes(self, absolute=False):
    249235        """Return all nodes in mesh.
    … …  
    256242        are to be made absolute by taking georeference into account
    257243        Default is False as many parts of ANUGA expects relative coordinates.
    258         (To see which, switch to default absolute=True and run tests).       
     244        (To see which, switch to default absolute=True and run tests).
    259245        """
    260246
    … …  
    264250            if not self.geo_reference.is_absolute():
    265251                V = self.geo_reference.get_absolute(V)
    266                
     252
    267253        return V
    268    
    269     def get_node(self, i,
    270                  absolute=False):
     254
     255    def get_node(self, i, absolute=False):
    271256        """Return node coordinates for triangle i.
    272257
    … …  
    274259        are to be made absolute by taking georeference into account
    275260        Default is False as many parts of ANUGA expects relative coordinates.
    276         (To see which, switch to default absolute=True and run tests).       
    277         """
    278 
    279        
     261        (To see which, switch to default absolute=True and run tests).
     262
     263        Note: This method returns a modified _copy_ of the nodes slice if
     264              absolute is True.  If absolute is False, just return the slice.
     265              This is related to the ensure_numeric() returning a copy problem.
     266        """
     267
    280268        V = self.nodes[i,:]
    281269        if absolute is True:
    282270            if not self.geo_reference.is_absolute():
    283                 return V + num.array([self.geo_reference.get_xllcorner(),
    284                                       self.geo_reference.get_yllcorner()], num.Float)
    285             else:
    286                 return V
    287         else:
    288             return V
    289    
    290        
    291 
    292     def get_vertex_coordinates(self,
    293                                triangle_id=None,
    294                                absolute=False):
    295         """Return vertex coordinates for all triangles.
    296        
     271                # get a copy so as not to modify the internal self.nodes array
     272                V = copy.copy(V)
     273                V += num.array([self.geo_reference.get_xllcorner(),
     274                                self.geo_reference.get_yllcorner()], num.float)
     275        return V
     276
     277    def get_vertex_coordinates(self, triangle_id=None, absolute=False):
     278        """Return vertex coordinates for all triangles.
     279
    297280        Return all vertex coordinates for all triangles as a 3*M x 2 array
    298281        where the jth vertex of the ith triangle is located in row 3*i+j and
    … …  
    301284        if triangle_id is specified (an integer) the 3 vertex coordinates
    302285        for triangle_id are returned.
    303        
     286
    304287        Boolean keyword argument absolute determines whether coordinates
    305288        are to be made absolute by taking georeference into account
    306289        Default is False as many parts of ANUGA expects relative coordinates.
    307290        """
    308        
     291
    309292        V = self.vertex_coordinates
    310293
    311         if triangle_id is None:   
     294        if triangle_id is None:
    312295            if absolute is True:
    313296                if not self.geo_reference.is_absolute():
    314297                    V = self.geo_reference.get_absolute(V)
    315      
    316298            return V
    317299        else:
    … …  
    320302            assert int(i) == i, msg
    321303            assert 0 <= i < self.number_of_triangles
    322            
    323             i3 = 3*i 
     304
     305            i3 = 3*i
    324306            if absolute is True and not self.geo_reference.is_absolute():
    325307                offset=num.array([self.geo_reference.get_xllcorner(),
    326                                   self.geo_reference.get_yllcorner()], num.Float)
     308                                  self.geo_reference.get_yllcorner()], num.float)
    327309                return num.array([V[i3,:]+offset,
    328310                                  V[i3+1,:]+offset,
    329                                   V[i3+2,:]+offset], num.Float)
     311                                  V[i3+2,:]+offset], num.float)
    330312            else:
    331                 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.Float)
    332                
    333 
     313                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
    334314
    335315    def get_vertex_coordinate(self, i, j, absolute=False):
    … …  
    340320        msg = 'vertex id j must be an integer in [0,1,2]'
    341321        assert j in [0,1,2], msg
    342        
    343         V = self.get_vertex_coordinates(triangle_id=i,
    344                                         absolute=absolute)
     322
     323        V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute)
    345324        return V[j,:]
    346    
    347 
    348325
    349326    def compute_vertex_coordinates(self):
    … …  
    356333
    357334        M = self.number_of_triangles
    358         vertex_coordinates = num.zeros((3*M, 2), num.Float)
     335        vertex_coordinates = num.zeros((3*M, 2), num.float)
    359336
    360337        for i in range(M):
    … …  
    365342        return vertex_coordinates
    366343
    367 
    368 
    369344    def get_triangles(self, indices=None):
    370345        """Get mesh triangles.
    … …  
    382357        if indices is None:
    383358            return self.triangles
    384             #indices = range(M)
    385 
    386         return num.take(self.triangles, indices)
    387    
    388 
     359
     360        return num.take(self.triangles, indices, axis=0)
    389361
    390362    def get_disconnected_triangles(self):
    … …  
    405377
    406378        The triangles created will have the format
    407 
    408379        [[0,1,2],
    409380         [3,4,5],
    410381         [6,7,8],
    411382         ...
    412          [3*M-3 3*M-2 3*M-1]]         
     383         [3*M-3 3*M-2 3*M-1]]
    413384        """
    414385
    415386        M = len(self) # Number of triangles
    416387        K = 3*M       # Total number of unique vertices
    417        
    418         T = num.reshape(num.arange(K, typecode=num.Int), (M,3))
    419        
    420         return T     
    421 
    422    
    423 
    424     def get_unique_vertices(self,  indices=None):
    425         """FIXME(Ole): This function needs a docstring
    426         """
     388        return num.reshape(num.arange(K, dtype=num.int), (M,3))
     389
     390    def get_unique_vertices(self, indices=None):
     391        """FIXME(Ole): This function needs a docstring"""
     392
    427393        triangles = self.get_triangles(indices=indices)
    428394        unique_verts = {}
    429395        for triangle in triangles:
    430            unique_verts[triangle[0]] = 0
    431            unique_verts[triangle[1]] = 0
    432            unique_verts[triangle[2]] = 0
     396            unique_verts[triangle[0]] = 0
     397            unique_verts[triangle[1]] = 0
     398            unique_verts[triangle[2]] = 0
    433399        return unique_verts.keys()
    434 
    435400
    436401    def get_triangles_and_vertices_per_node(self, node=None):
    … …  
    447412            # Get index for this node
    448413            first = num.sum(self.number_of_triangles_per_node[:node])
    449            
     414
    450415            # Get number of triangles for this node
    451416            count = self.number_of_triangles_per_node[node]
    … …  
    459424                triangle_list.append( (volume_id, vertex_id) )
    460425
    461             triangle_list = num.array(triangle_list, num.Int)    #array default#
     426            triangle_list = num.array(triangle_list, num.int)    #array default#
    462427        else:
    463428            # Get info for all nodes recursively.
    … …  
    465430            # working directly with the inverted triangle structure
    466431            for i in range(self.number_of_full_nodes):
    467                
    468432                L = self.get_triangles_and_vertices_per_node(node=i)
    469 
    470433                triangle_list.append(L)
    471434
    472435        return triangle_list
    473        
    474436
    475437    def build_inverted_triangle_structure(self):
    … …  
    481443        listing for each node how many triangles use it. N is the number of
    482444        nodes in mesh.
    483        
     445
    484446        vertex_value_indices: An array of length M listing indices into
    485447        triangles ordered by node number. The (triangle_id, vertex_id)
    … …  
    489451        used to average vertex values efficiently.
    490452
    491        
    492453        Example:
    493        
    494454        a = [0.0, 0.0] # node 0
    495455        b = [0.0, 2.0] # node 1
    … …  
    498458        e = [2.0, 2.0] # node 4
    499459        f = [4.0, 0.0] # node 5
    500 
    501460        nodes = array([a, b, c, d, e, f])
    502        
    503         #bac, bce, ecf, dbe, daf, dae
    504         triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])       
    505 
    506 
    507         For this structure
    508 
     461
     462        #                    bac,     bce,     ecf,     dbe
     463        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
     464
     465        For this structure:
    509466        number_of_triangles_per_node = [1 3 3 1 3 1]
    510467        which means that node a has 1 triangle associated with it, node b
    511468        has 3, node has 3 and so on.
    512        
     469
    513470        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
    514471        which reflects the fact that
    … …  
    524481                   and by triangle 2, vertex 0 (index = 6)
    525482                   and by triangle 3, vertex 2 (index = 11)
    526         node 5 is used by triangle 2, vertex 2 (index = 8)                   
    527        
     483        node 5 is used by triangle 2, vertex 2 (index = 8)
    528484
    529485        Preconditions:
    … …  
    532488        Postcondition:
    533489          self.number_of_triangles_per_node is built
    534           self.vertex_value_indices is built         
     490          self.vertex_value_indices is built
    535491        """
    536492
    537493        # Count number of triangles per node
    538         number_of_triangles_per_node = num.zeros(self.number_of_full_nodes)
     494        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes,
     495                                                 num.int)       #array default#
    539496        for volume_id, triangle in enumerate(self.get_triangles()):
    540497            for vertex_id in triangle:
    … …  
    543500        # Allocate space for inverted structure
    544501        number_of_entries = num.sum(number_of_triangles_per_node)
    545         vertex_value_indices = num.zeros(number_of_entries)
     502        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
    546503
    547504        # Register (triangle, vertex) indices for each node
    548         vertexlist = [None]*self.number_of_full_nodes
     505        vertexlist = [None] * self.number_of_full_nodes
    549506        for volume_id in range(self.number_of_full_triangles):
    550 
    551507            a = self.triangles[volume_id, 0]
    552508            b = self.triangles[volume_id, 1]
    553509            c = self.triangles[volume_id, 2]
    554510
    555             for vertex_id, node_id in enumerate([a,b,c]):
     511            for vertex_id, node_id in enumerate([a, b, c]):
    556512                if vertexlist[node_id] is None:
    557513                    vertexlist[node_id] = []
    558        
    559                 vertexlist[node_id].append( (volume_id, vertex_id) )
    560 
     514                vertexlist[node_id].append((volume_id, vertex_id))
    561515
    562516        # Build inverted triangle index array
    … …  
    566520                for volume_id, vertex_id in vertices:
    567521                    vertex_value_indices[k] = 3*volume_id + vertex_id
    568                                        
    569522                    k += 1
    570523
    … …  
    573526        self.vertex_value_indices = vertex_value_indices
    574527
    575 
    576        
    577 
    578528    def get_extent(self, absolute=False):
    579529        """Return min and max of all x and y coordinates
    … …  
    582532        are to be made absolute by taking georeference into account
    583533        """
    584        
    585 
    586534
    587535        C = self.get_vertex_coordinates(absolute=absolute)
    … …  
    589537        Y = C[:,1:6:2].copy()
    590538
    591         xmin = min(X.flat)
    592         xmax = max(X.flat)
    593         ymin = min(Y.flat)
    594         ymax = max(Y.flat)
    595         #print "C",C
     539        xmin = num.min(X)
     540        xmax = num.max(X)
     541        ymin = num.min(Y)
     542        ymax = num.max(Y)
     543
    596544        return xmin, xmax, ymin, ymax
    597545
    598546    def get_areas(self):
    599         """Get areas of all individual triangles.
    600         """
    601         return self.areas       
     547        """Get areas of all individual triangles."""
     548
     549        return self.areas
    602550
    603551    def get_area(self):
    604         """Return total area of mesh
    605         """
     552        """Return total area of mesh"""
    606553
    607554        return num.sum(self.areas)
    … …  
    609556    def set_georeference(self, g):
    610557        self.geo_reference = g
    611        
     558
    612559    def get_georeference(self):
    613560        return self.geo_reference
    614        
    615        
    616                
     561

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

    r6928 r7276  
    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):
    … …  
    129129
    130130        try:
    131             q = num.array(q, num.Float)
     131            q = num.array(q, num.float)
    132132        except:
    133133            msg = 'Return value from time boundary function could '
    134             msg += 'not be converted into a Numeric array of floats.\n'
     134            msg += 'not be converted into a numeric array of floats.\n'
    135135            msg += 'Specified function should return either list or array.\n'
    136136            msg += 'I got %s' %str(q)
    … …  
    241241
    242242        if verbose: print 'Find midpoint coordinates of entire boundary'
    243         self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.Float)
     243        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
    244244        boundary_keys = domain.boundary.keys()
    245245
    … …  
    261261           
    262262            # Compute midpoints
    263             if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
    264             if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
    265             if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
     263            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
     264            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
     265            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
    266266
    267267            # Convert to absolute UTM coordinates
    … …  
    371371                        # See http://docs.python.org/lib/warning-filter.html
    372372                        self.default_boundary_invoked = True
    373                    
    374 
    375             if res == NAN:
     373           
     374            if num.any(res == NAN):
    376375                x,y=self.midpoint_coordinates[i,:]
    377376                msg = 'NAN value found in file_boundary at '
    … …  
    423422
    424423    This was added by Nils Goseberg et al in April 2009
    425        
    426424    """
    427425
    … …  
    429427                 use_cache=False, verbose=False):
    430428        import time
    431         from Numeric import array, zeros, Float
    432429        from anuga.config import time_format
    433430        from anuga.abstract_2d_finite_volumes.util import file_function
    … …  
    445442
    446443        if verbose: print 'Find midpoint coordinates of entire boundary'
    447         self.midpoint_coordinates = zeros( (len(domain.boundary), 2), Float)
     444        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
    448445        boundary_keys = domain.boundary.keys()
    449 
    450446
    451447        xllcorner = domain.geo_reference.get_xllcorner()
    452448        yllcorner = domain.geo_reference.get_yllcorner()       
    453        
    454449
    455450        # Make ordering unique #FIXME: should this happen in domain.py?
    … …  
    459454        self.boundary_indices = {}
    460455        for i, (vol_id, edge_id) in enumerate(boundary_keys):
    461 
    462456            base_index = 3*vol_id
    463457            x0, y0 = V[base_index, :]
    … …  
    466460           
    467461            # Compute midpoints
    468             if edge_id == 0: m = array([(x1 + x2)/2, (y1 + y2)/2])
    469             if edge_id == 1: m = array([(x0 + x2)/2, (y0 + y2)/2])
    470             if edge_id == 2: m = array([(x1 + x0)/2, (y1 + y0)/2])
     462            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2])
     463            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2])
     464            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2])
    471465
    472466            # Convert to absolute UTM coordinates
    … …  
    483477        if verbose: print 'Initialise file_function'
    484478        self.F = file_function(filename, domain,
    485                                interpolation_points=self.midpoint_coordinates,
     479                                   interpolation_points=self.midpoint_coordinates,
    486480                               time_thinning=time_thinning,
    487481                               use_cache=use_cache,
    … …  
    497491        # only part of mesh boundary is actually used with a given
    498492        # file boundary sww file.
    499         if hasattr(self.F, 'indices_outside_mesh') and\
    500                len(self.F.indices_outside_mesh) > 0:
     493        if (hasattr(self.F, 'indices_outside_mesh') and
     494               len(self.F.indices_outside_mesh)) > 0:
    501495            msg = 'WARNING: File_boundary has points outside the mesh '
    502             msg += 'given in %s. ' %filename
     496            msg += 'given in %s. ' % filename
    503497            msg += 'See warning message issued by Interpolation_function '
    504498            msg += 'for details (should appear above somewhere if '
    … …  
    510504            #raise Exception(msg)
    511505
    512        
    513506        q = self.F(0, point_id=0)
    514507
    515508        d = len(domain.conserved_quantities)
    516         msg = 'Values specified in file %s must be ' %filename
    517         msg += ' a list or an array of length %d' %d
     509        msg = 'Values specified in file %s must be ' % filename
     510        msg += 'a list or an array of length %d' % d
    518511        assert len(q) == d, msg
    519512
    … …  
    525518    def evaluate(self, vol_id=None, edge_id=None):
    526519        """Return linearly interpolated values based on domain.time
    527         at midpoint of segment defined by vol_id and edge_id.
     520               at midpoint of segment defined by vol_id and edge_id.
    528521        """
    529         q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
     522
     523        q = self.domain.get_conserved_quantities(vol_id, edge=edge_id)
    530524        t = self.domain.time
    531525
    532526        if vol_id is not None and edge_id is not None:
    533             i = self.boundary_indices[ vol_id, edge_id ]
    534             res = self.F(t, point_id = i)
     527            i = self.boundary_indices[vol_id, edge_id]
     528            res = self.F(t, point_id=i)
    535529
    536530            if res == NAN:
    537                 x,y=self.midpoint_coordinates[i,:]
     531                x,y = self.midpoint_coordinates[i,:]
    538532                msg = 'NAN value found in file_boundary at '
    539                 msg += 'point id #%d: (%.2f, %.2f).\n' %(i, x, y)
    540 
    541                 if hasattr(self.F, 'indices_outside_mesh') and\
    542                        len(self.F.indices_outside_mesh) > 0:
     533                msg += 'point id #%d: (%.2f, %.2f).\n' % (i, x, y)
     534
     535                if (hasattr(self.F, 'indices_outside_mesh') and
     536                       len(self.F.indices_outside_mesh) > 0):
    543537                    # Check if NAN point is due it being outside
    544538                    # boundary defined in sww file.
    … …  
    546540                    if i in self.F.indices_outside_mesh:
    547541                        msg += 'This point refers to one outside the '
    548                         msg += 'mesh defined by the file %s.\n'\
    549                                %self.F.filename
     542                        msg += 'mesh defined by the file %s.\n' % self.F.filename
    550543                        msg += 'Make sure that the file covers '
    551544                        msg += 'the boundary segment it is assigned to '
    … …  
    553546                    else:
    554547                        msg += 'This point is inside the mesh defined '
    555                         msg += 'the file %s.\n' %self.F.filename
     548                        msg += 'the file %s.\n' % self.F.filename
    556549                        msg += 'Check this file for NANs.'
    557550                raise Exception, msg

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

    r6717 r7276  
    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
    … …  
    203203
    204204
    205 
    206205def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
    207206
    … …  
    260259    E = EIndex(n,m)
    261260
    262     points = num.zeros( (Np,2), num.Float)
     261    points = num.zeros( (Np,2), num.float)
    263262
    264263    for i in range(m+1):
    … …  
    272271
    273272
    274     elements = num.zeros( (Nt,3), num.Int)
     273    elements = num.zeros( (Nt,3), num.int)
    275274    boundary = {}
    276275    Idgl = []
    … …  
    340339    Idgr.extend(Idgl)
    341340
    342     Idfl = num.array(Idfl, num.Int)
    343     Idgr = num.array(Idgr, num.Int)
     341    Idfl = num.array(Idfl, num.int)
     342    Idgr = num.array(Idgr, num.int)
    344343   
    345344    full_send_dict[processor]  = [Idfl, Idfl]

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

    r6654 r7276  
    1212from math import pi, sqrt
    1313
    14 import Numeric as num
    15        
     14import numpy as num
     15
    1616
    1717class Mesh(General_mesh):
    … …  
    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
    … …  
    9191                              verbose=verbose)
    9292
    93         if verbose: print 'Initialising mesh'         
     93        if verbose: print 'Initialising mesh'
    9494
    9595        N = len(self) #Number_of_triangles
    … …  
    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
    … …  
    111111
    112112        #Initialise each triangle
    113         if verbose: print 'Mesh: Computing centroids and radii'       
     113        if verbose: print 'Mesh: Computing centroids and radii'
    114114        for i in range(N):
    115115            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
    … …  
    117117            x0, y0 = V[3*i, :]
    118118            x1, y1 = V[3*i+1, :]
    119             x2, y2 = V[3*i+2, :]                       
     119            x2, y2 = V[3*i+2, :]
    120120
    121121            #x0 = V[i, 0]; y0 = V[i, 1]
    … …  
    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
    … …  
    166166
    167167        #Build neighbour structure
    168         if verbose: print 'Mesh: Building neigbour structure'               
     168        if verbose: print 'Mesh: Building neigbour structure'
    169169        self.build_neighbour_structure()
    170170
    … …  
    178178
    179179        #Build tagged element  dictionary mapping (tag) to array of elements
    180         if verbose: print 'Mesh: Building tagged elements dictionary'       
     180        if verbose: print 'Mesh: Building tagged elements dictionary'
    181181        self.build_tagged_elements_dictionary(tagged_elements)
    182182
    … …  
    184184        self.lone_vertices = []
    185185        #Check that all vertices have been registered
    186         for node, count in enumerate(self.number_of_triangles_per_node):   
     186        for node, count in enumerate(self.number_of_triangles_per_node):
    187187            #msg = 'Node %d does not belong to an element.' %node
    188188            #assert count > 0, msg
    189189            if count == 0:
    190190                self.lone_vertices.append(node)
    191                
     191
    192192        #Update boundary indices FIXME: OBSOLETE
    193193        #self.build_boundary_structure()
    194194
    195195        #FIXME check integrity?
    196         if verbose: print 'Mesh: Done'               
     196        if verbose: print 'Mesh: Done'
    197197
    198198    def __repr__(self):
    … …  
    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. '
    … …  
    400400        #print "tagged_elements", tagged_elements
    401401        self.tagged_elements = tagged_elements
    402        
     402
    403403    def get_tagged_elements(self):
    404404        return self.tagged_elements
    … …  
    459459
    460460        Using the mesh boundary, derive a bounding polygon for this mesh.
    461         If multiple vertex values are present (vertices stored uniquely), 
    462         the algorithm will select the path that contains the entire mesh.
     461        If multiple vertex values are present (vertices stored uniquely),
     462        the algorithm will select the path that contains the entire mesh.
    463463
    464464        All points are in absolute UTM coordinates
    465465        """
    466        
    467         from anuga.utilities.numerical_tools import angle, ensure_numeric     
    468 
     466
     467        from anuga.utilities.numerical_tools import angle, ensure_numeric
    469468
    470469        # Get mesh extent
    471470        xmin, xmax, ymin, ymax = self.get_extent(absolute=True)
    472471        pmin = ensure_numeric([xmin, ymin])
    473         pmax = ensure_numeric([xmax, ymax])       
    474 
     472        pmax = ensure_numeric([xmax, ymax])
    475473
    476474        # Assemble dictionary of boundary segments and choose starting point
    … …  
    478476        inverse_segments = {}
    479477        p0 = None
    480         mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh
     478
     479        # Start value across entire mesh
     480        mindist = num.sqrt(num.sum((pmax-pmin)**2))
    481481        for i, edge_id in self.boundary.keys():
    482482            # Find vertex ids for boundary segment
    … …  
    485485            if edge_id == 2: a = 0; b = 1
    486486
    487             A = self.get_vertex_coordinate(i, a, absolute=True) # Start
    488             B = self.get_vertex_coordinate(i, b, absolute=True) # End
    489 
     487            A = self.get_vertex_coordinate(i, a, absolute=True)    # Start
     488            B = self.get_vertex_coordinate(i, b, absolute=True)    # End
    490489
    491490            # Take the point closest to pmin as starting point
    … …  
    503502                p0 = B
    504503
    505 
    506504            # Sanity check
    507505            if p0 is None:
    508                 raise Exception('Impossible')
    509 
     506                msg = 'Impossible: p0 is None!?'
     507                raise Exception, msg
    510508
    511509            # Register potential paths from A to B
    512510            if not segments.has_key(tuple(A)):
    513                 segments[tuple(A)] = [] # Empty list for candidate points
    514 
    515             segments[tuple(A)].append(B)               
    516 
     511                segments[tuple(A)] = []    # Empty list for candidate points
     512
     513            segments[tuple(A)].append(B)
    517514
    518515        # Start with smallest point and follow boundary (counter clock wise)
    519516        polygon = [list(p0)]# Storage for final boundary polygon
    520         point_registry = {} # Keep track of storage to avoid multiple runs 
    521                             # around boundary. This will only be the case if
    522                             # there are more than one candidate.
     517        point_registry = {} # Keep track of storage to avoid multiple runs
     518                            # around boundary. This will only be the case if
     519                            # there are more than one candidate.
    523520                            # FIXME (Ole): Perhaps we can do away with polygon
    524521                            # and use only point_registry to save space.
    525522
    526         point_registry[tuple(p0)] = 0                           
    527                            
     523        point_registry[tuple(p0)] = 0
     524
    528525        while len(point_registry) < len(self.boundary):
    529 
    530526            candidate_list = segments[tuple(p0)]
    531527            if len(candidate_list) > 1:
    532                 # Multiple points detected (this will be the case for meshes 
    533                 # with duplicate points as those used for discontinuous
    534                 # triangles with vertices stored uniquely).
    535                 # Take the candidate that is furthest to the clockwise 
    536                 # direction, as that will follow the boundary.
    537                 #
    538                 # This will also be the case for pathological triangles
    539                 # that have no neighbours.
     528                # Multiple points detected (this will be the case for meshes
     529                # with duplicate points as those used for discontinuous
     530                # triangles with vertices stored uniquely).
     531                # Take the candidate that is furthest to the clockwise
     532                # direction, as that will follow the boundary.
     533                #
     534                # This will also be the case for pathological triangles
     535                # that have no neighbours.
    540536
    541537                if verbose:
    542                     print 'Point %s has multiple candidates: %s'\
    543                           %(str(p0), candidate_list)
     538                    print ('Point %s has multiple candidates: %s'
     539                           % (str(p0), candidate_list))
    544540
    545541                # Check that previous are not in candidate list
    … …  
    548544
    549545                # Choose vector against which all angles will be measured
    550                 if len(polygon) > 1:   
    551                     v_prev = p0 - polygon[-2] # Vector that leads to p0
    552                                               # from previous point
     546                if len(polygon) > 1:
     547                    v_prev = p0 - polygon[-2]    # Vector that leads to p0
     548                                                 # from previous point
    553549                else:
    554                     # FIXME (Ole): What do we do if the first point has 
    555                     # multiple candidates?
     550                    # FIXME (Ole): What do we do if the first point has
     551                    # multiple candidates?
    556552                    # Being the lower left corner, perhaps we can use the
    557                     # vector [1, 0], but I really don't know if this is 
    558                     # completely watertight.
     553                    # vector [1, 0], but I really don't know if this is
     554                    # completely watertight.
    559555                    v_prev = [1.0, 0.0]
    560                    
    561 
    562                 # Choose candidate with minimum angle   
     556
     557                # Choose candidate with minimum angle
    563558                minimum_angle = 2*pi
    564559                for pc in candidate_list:
    565 
    566                     vc = pc-p0  # Candidate vector (from p0 to candidate pt)
    567                    
     560                    vc = pc-p0    # Candidate vector (from p0 to candidate pt)
     561
    568562                    # Angle between each candidate and the previous vector
    569563                    # in [-pi, pi]
    570564                    ac = angle(vc, v_prev)
    571565                    if ac > pi:
    572                         # Give preference to angles on the right hand side 
    573                         # of v_prev
     566                        # Give preference to angles on the right hand side
     567                        # of v_prev
    574568                        # print 'pc = %s, changing angle from %f to %f'\
    575                         # %(pc, ac*180/pi, (ac-2*pi)*180/pi)
     569                        # %(pc, ac*180/pi, (ac-2*pi)*180/pi)
    576570                        ac = ac-2*pi
    577571
    578                     # Take the minimal angle corresponding to the 
    579                     # rightmost vector
     572                    # Take the minimal angle corresponding to the
     573                    # rightmost vector
    580574                    if ac < minimum_angle:
    581575                        minimum_angle = ac
    582                         p1 = pc             # Best candidate
    583                        
     576                        p1 = pc             # Best candidate
    584577
    585578                if verbose is True:
    586579                    print '  Best candidate %s, angle %f'\
    587580                          %(p1, minimum_angle*180/pi)
    588                    
    589581            else:
    590582                p1 = candidate_list[0]
    591583
    592                
    593584            if point_registry.has_key(tuple(p1)):
    594                 # We have reached a point already visited.
    595                
    596                 if num.allclose(p1, polygon[0]):
    597                     # If it is the initial point, the polygon is complete.
    598                    
     585                # We have reached a point already visited.
     586                if num.allclose(p1, polygon[0]):
     587                    # If it is the initial point, the polygon is complete.
    599588                    if verbose is True:
    600                         print '  Stop criterion fulfilled at point %s' %p1
    601                         print polygon               
    602                        
     589                        print '  Stop criterion fulfilled at point %s' %p1
     590                        print polygon
     591
    603592                    # We have completed the boundary polygon - yeehaa
    604                     break
    605                 else:   
    606                     # The point already visited is not the initial point
    607                     # This would be a pathological triangle, but the
    608                     # algorithm must be able to deal with this
    609                     pass
    610    
     593                    break
     594                else:
     595                    # The point already visited is not the initial point
     596                    # This would be a pathological triangle, but the
     597                    # algorithm must be able to deal with this
     598                    pass
     599
    611600            else:
    612                 # We are still finding new points on the boundary
     601                # We are still finding new points on the boundary
    613602                point_registry[tuple(p1)] = len(point_registry)
    614            
    615             polygon.append(list(p1)) # De-Numeric each point :-)
     603
     604            polygon.append(list(p1))    # De-numeric each point :-)
    616605            p0 = p1
    617606
    618 
    619607        return polygon
    620 
    621608
    622609    def check_integrity(self):
    … …  
    641628            x1, y1 = V[3*i+1, :]
    642629            x2, y2 = V[3*i+2, :]
    643            
     630
    644631            # Check that area hasn't been compromised
    645632            area = self.areas[i]
    … …  
    678665
    679666                x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product
    680                
     667
    681668                msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v)
    682669                msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,))
    683                 msg += ' Inner product is %e.' %x               
     670                msg += ' Inner product is %e.' %x
    684671                assert x < epsilon, msg
    685672
    … …  
    690677        for i in range(N):
    691678            # For each triangle
    692            
     679
    693680            for k, neighbour_id in enumerate(self.neighbours[i,:]):
    694681
    … …  
    752739                # Node is lone - i.e. not part of the mesh
    753740                continue
    754            
     741
    755742            k += 1
    756            
     743
    757744            volume_id = index / 3
    758745            vertex_id = index % 3
    759            
     746
    760747            msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\
    761748                  %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node)
    762749            assert self.triangles[volume_id, vertex_id] == current_node, msg
    763                        
     750
    764751            if self.number_of_triangles_per_node[current_node] == k:
    765752                # Move on to next node
    … …  
    788775            if not self.geo_reference.is_absolute():
    789776                V = self.geo_reference.get_absolute(V)
    790            
     777
    791778        return V
    792779
    793        
     780
    794781    def get_radii(self):
    795782        """Return all radii.
    796783        Return radius of inscribed cirle for all triangles
    797784        """
    798         return self.radii       
     785        return self.radii
    799786
    800787
    … …  
    826813        str += '  Areas [m^2]:\n'
    827814        str += '    A in [%f, %f]\n' %(min(areas), max(areas))
    828         str += '    number of distinct areas: %d\n' %(len(areas))       
     815        str += '    number of distinct areas: %d\n' %(len(areas))
    829816        str += '    Histogram:\n'
    830817
    … …  
    833820            lo = hi
    834821            if i+1 < len(bins):
    835                 #Open upper interval               
     822                #Open upper interval
    836823                hi = bins[i+1]
    837                 str += '      [%f, %f[: %d\n' %(lo, hi, count)               
     824                str += '      [%f, %f[: %d\n' %(lo, hi, count)
    838825            else:
    839826                #Closed upper interval
    … …  
    849836            k = 0
    850837            lower = min(areas)
    851             for i, a in enumerate(areas):       
    852                 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas               
     838            for i, a in enumerate(areas):
     839                if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas
    853840                    str += '      %d triangles in [%f, %f]\n' %(i-k, lower, a)
    854841                    lower = a
    855842                    k = i
    856                    
     843
    857844            str += '      %d triangles in [%f, %f]\n'\
    858                    %(N-k, lower, max(areas))                   
    859                
    860                      
     845                   %(N-k, lower, max(areas))
     846
     847
    861848        str += '  Boundary:\n'
    862849        str += '    Number of boundary segments == %d\n' %(len(self.boundary))
    863         str += '    Boundary tags == %s\n' %self.get_boundary_tags() 
     850        str += '    Boundary tags == %s\n' %self.get_boundary_tags()
    864851        str += '------------------------------------------------\n'
    865        
     852
    866853
    867854        return str
    868    
     855
    869856
    870857    def get_triangle_containing_point(self, point):
    … …  
    874861
    875862        """
    876        
     863
    877864        # FIXME(Ole): This function is currently brute force
    878865        # because I needed it for diagnostics.
    … …  
    884871        polygon = self.get_boundary_polygon()
    885872        #print polygon, point
    886        
     873
    887874        if is_outside_polygon(point, polygon):
    888875            msg = 'Point %s is outside mesh' %str(point)
    … …  
    899886            if is_inside_polygon(point, poly, closed=True):
    900887                return i
    901                
     888
    902889        return
    903890
    … …  
    930917        Resulting segments are unsorted
    931918        """
    932        
     919
    933920        V = self.get_vertex_coordinates()
    934921        N = len(self)
    935        
     922
    936923        # Adjust polyline to mesh spatial origin
    937924        polyline = self.geo_reference.get_relative(polyline)
    … …  
    939926        if use_cache is True:
    940927            segments = cache(get_intersecting_segments,
    941                              (V, N, polyline), 
     928                             (V, N, polyline),
    942929                             {'verbose': verbose},
    943930                             verbose=verbose)
    944         else:                 
     931        else:
    945932            segments = get_intersecting_segments(V, N, polyline,
    946933                                                 verbose=verbose)
    947        
     934
    948935
    949936        return segments
    950        
    951  
     937
     938
    952939
    953940    def get_triangle_neighbours(self, tri_id):
    … …  
    956943
    957944        Negative returned triangle id's represent a boundary as a neighbour.
    958        
     945
    959946        If the given triangle id is bad, return an empty list.
    960947        """
    … …  
    964951        except IndexError:
    965952            return []
    966        
     953
    967954
    968955    def get_interpolation_object(self):
    969956        """Get object I that will allow linear interpolation using this mesh
    970        
    971         This is a time consuming process but it needs only to be 
     957
     958        This is a time consuming process but it needs only to be
    972959        once for the mesh.
    973        
    974         Interpolation can then be done using 
    975        
    976         result = I.interpolate_block(vertex_values, interpolation_points)       
    977        
     960
     961        Interpolation can then be done using
     962
     963        result = I.interpolate_block(vertex_values, interpolation_points)
     964
    978965        where vertex values have been obtained from a quantity using
    979966        vertex_values, triangles = self.get_vertex_values()
    … …  
    984971        else:
    985972            from anuga.fit_interpolate.interpolate import Interpolate
    986                        
    987             # Get discontinuous mesh - this will match internal 
     973
     974            # Get discontinuous mesh - this will match internal
    988975            # representation of vertex values
    989976            triangles = self.get_disconnected_triangles()
    … …  
    992979            I = Interpolate(vertex_coordinates, triangles)
    993980            self.interpolation_object = I
    994        
    995         return I               
    996        
     981
     982        return I
     983
    997984
    998985class Triangle_intersection:
    999986    """Store information about line segments intersecting a triangle
    1000    
     987
    1001988    Attributes are
    1002989
    … …  
    10141001                 length=None,
    10151002                 triangle_id=None):
    1016         self.segment = segment             
     1003        self.segment = segment
    10171004        self.normal = normal
    10181005        self.length = length
    10191006        self.triangle_id = triangle_id
    1020        
     1007
    10211008
    10221009    def __repr__(self):
    … …  
    10271014               self.length,
    10281015               self.triangle_id)
    1029    
     1016
    10301017        return s
    10311018
    1032        
     1019
    10331020
    10341021def _get_intersecting_segments(V, N, line,
    … …  
    10511038    from anuga.utilities.polygon import intersection
    10521039    from anuga.utilities.polygon import is_inside_polygon
    1053  
     1040
    10541041    msg = 'Line segment must contain exactly two points'
    10551042    assert len(line) == 2, msg
    … …  
    10601047    eta0 = line[0][1]
    10611048
    1062  
     1049
    10631050    # Check intersection with edge segments for all triangles
    10641051    # FIXME (Ole): This should be implemented in C
    … …  
    10691056        x1, y1 = V[3*i+1, :]
    10701057        x2, y2 = V[3*i+2, :]
    1071          
     1058
    10721059
    10731060        edge_segments = [[[x0,y0], [x1, y1]],
    … …  
    10761063
    10771064        # Find segments that are intersected by line
    1078        
     1065
    10791066        intersections = {} # Use dictionary to record points only once
    10801067        for edge in edge_segments:
    … …  
    10821069            status, value = intersection(line, edge)
    10831070            #if value is not None: print 'Triangle %d, Got' %i, status, value
    1084                
     1071
    10851072            if status == 1:
    10861073                # Normal intersection of one edge or vertex
    1087                 intersections[tuple(value)] = i                 
     1074                intersections[tuple(value)] = i
    10881075
    10891076                # Exclude singular intersections with vertices
    … …  
    11011088                # along this edge, it will be recorded here.
    11021089                intersections[tuple(value[0,:])] = i
    1103                 intersections[tuple(value[1,:])] = i                                   
    1104 
    1105                
     1090                intersections[tuple(value[1,:])] = i
     1091
     1092
    11061093        if len(intersections) == 1:
    11071094            # Check if either line end point lies fully within this triangle
    … …  
    11131100                intersections[tuple(line[1])] = i
    11141101            elif is_inside_polygon(line[0], poly, closed=False):
    1115                 intersections[tuple(line[0])] = i         
     1102                intersections[tuple(line[0])] = i
    11161103            else:
    1117                 # Ignore situations where one vertex is touch, for instance                   
     1104                # Ignore situations where one vertex is touch, for instance
    11181105                continue
    11191106
    … …  
    11391126
    11401127            # 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)
    1143             d0 = num.sqrt(num.sum(z0**2)) 
     1128            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
     1129            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
     1130            d0 = num.sqrt(num.sum(z0**2))
    11441131            d1 = num.sqrt(num.sum(z1**2))
    1145                
     1132
    11461133            if d1 < d0:
    11471134                # Swap
    … …  
    11511138
    11521139            # (x0,y0) is now the origin of the intersecting segment
    1153                
     1140
    11541141
    11551142            # Normal direction:
    11561143            # Right hand side relative to line direction
    1157             vector = num.array([x1 - x0, y1 - y0], num.Float) # Segment vector
     1144            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
    11581145            length = num.sqrt(num.sum(vector**2))      # Segment length
    1159             normal = num.array([vector[1], -vector[0]], num.Float)/length
    1160 
    1161 
    1162             segment = ((x0,y0), (x1, y1))   
     1146            normal = num.array([vector[1], -vector[0]], num.float)/length
     1147
     1148
     1149            segment = ((x0,y0), (x1, y1))
    11631150            T = Triangle_intersection(segment=segment,
    11641151                                      normal=normal,
    … …  
    11681155
    11691156            # Add segment unless it was done earlier
    1170             if not triangle_intersections.has_key(segment):   
     1157            if not triangle_intersections.has_key(segment):
    11711158                triangle_intersections[segment] = T
    11721159
    11731160
    1174     # Return segments as a list           
     1161    # Return segments as a list
    11751162    return triangle_intersections.values()
    11761163
    11771164
    11781165def get_intersecting_segments(V, N, polyline,
    1179                               verbose=False):       
     1166                              verbose=False):
    11801167    """Internal function to find edges intersected by Polyline
    1181    
     1168
    11821169    Input:
    11831170        V: Vertex coordinates as obtained by mesh.get_vertex_coordinates()
    11841171        N: Number of triangles in mesh
    1185         polyline - list of points forming a segmented line       
     1172        polyline - list of points forming a segmented line
    11861173        verbose
    11871174    Output:
    … …  
    11901177    This method is used by the public method
    11911178    get_intersecting_segments(self, polyline) which also contains
    1192     more documentation.   
     1179    more documentation.
    11931180    """
    11941181
    11951182    msg = 'Polyline must contain at least two points'
    11961183    assert len(polyline) >= 2, msg
    1197      
    1198      
     1184
     1185
    11991186    # For all segments in polyline
    12001187    triangle_intersections = []
    … …  
    12061193            print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0[0], point0[1],
    12071194                                                  point1[0], point1[1])
    1208            
     1195
    12091196        line = [point0, point1]
    12101197        triangle_intersections += _get_intersecting_segments(V, N, line,
    … …  
    12141201    msg = 'No segments found'
    12151202    assert len(triangle_intersections) > 0, msg
    1216      
    1217        
     1203
     1204
    12181205    return triangle_intersections
    1219  
    1220 
    1221    
    1222            
    1223        
     1206
     1207
     1208
     1209
     1210
    12241211def segment_midpoints(segments):
    12251212    """Calculate midpoints of all segments
    1226    
     1213
    12271214    Inputs:
    12281215       segments: List of instances of class Segment
    1229        
     1216
    12301217    Ouputs:
    12311218       midpoints: List of points
    12321219    """
    1233    
     1220
    12341221    midpoints = []
    12351222    msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection'
    12361223    for segment in segments:
    12371224        assert isinstance(segment, Triangle_intersection), msg
    1238        
    1239         midpoint = num.sum(num.array(segment.segment, num.Float))/2
     1225
     1226        midpoint = num.sum(num.array(segment.segment, num.float), axis=0)/2
    12401227        midpoints.append(midpoint)
    12411228
    12421229    return midpoints
    1243    
    1244    
    1245    
     1230
     1231
     1232

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

    r6688 r7276  
    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
    161     that can be used to instanciate a domain object.
     166    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
     167    that can be used to instantiate a domain object.
    162168    """
    163169
    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 = {}
     182
     183    # num.transpose(None) gives scalar array of value None
    174184    point_atts = mesh_dict['vertex_attributes']
    175     point_titles  = mesh_dict['vertex_attribute_titles']
    176     geo_reference  = mesh_dict['geo_reference']
     185
     186    point_titles = mesh_dict['vertex_attribute_titles']
     187    geo_reference = mesh_dict['geo_reference']
    177188    if point_atts is not None:
    178         point_atts = num.transpose(point_atts)   
    179         for quantity, value_vector in map (None, point_titles, point_atts):
     189        point_atts = num.transpose(point_atts)
     190        for quantity, value_vector in map(None, point_titles, point_atts):
    180191            vertex_quantity_dict[quantity] = value_vector
    181192    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
    182193    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    183     return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
    184 
     194
     195    return (vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,
     196            tagged_elements_dict, geo_reference)
    185197
    186198
    187199def build_tagged_elements_dictionary(mesh_dict):
    188200    """Build the dictionary of element tags.
     201
    189202    tagged_elements is a dictionary of element arrays,
    190     keyed by tag:
    191     { (tag): [e1, e2, e3..] }
    192     """
     203    keyed by tag: { (tag): [e1, e2, e3..] }
     204    """
     205
    193206    tri_atts = mesh_dict['triangle_tags']
    194207    tagged_elements = {}
    … …  
    199212            tagged_elements.setdefault(tri_atts[tri_att_index],
    200213                                       []).append(tri_att_index)
     214
    201215    return tagged_elements
     216
    202217
    203218def pmesh_dict_to_tag_dict(mesh_dict):
    … …  
    205220    to a dictionary of tags, indexed with volume id and face number.
    206221    """
     222
    207223    triangles = mesh_dict['triangles']
    208224    sides = calc_sides(triangles)
    209225    tag_dict = {}
    210     for seg, tag in map(None, mesh_dict['segments'],
    211                         mesh_dict['segment_tags']):
     226    for seg, tag in map(None, mesh_dict['segments'], mesh_dict['segment_tags']):
    212227        v1 = int(seg[0])
    213228        v2 = int(seg[1])
    … …  
    223238
    224239def calc_sides(triangles):
    225     #Build dictionary mapping from sides (2-tuple of points)
    226     #to left hand side neighbouring triangle
     240    '''Build dictionary mapping from sides (2-tuple of points)
     241    to left hand side neighbouring triangle
     242    '''
     243
    227244    sides = {}
     245
    228246    for id, triangle in enumerate(triangles):
    229247        a = int(triangle[0])
    230248        b = int(triangle[1])
    231249        c = int(triangle[2])
     250
    232251        sides[a,b] = (id, 2) #(id, face)
    233252        sides[b,c] = (id, 0) #(id, face)
    234253        sides[c,a] = (id, 1) #(id, face)
     254
    235255    return sides
     256

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

    r7044 r7276  
    1414   Otherwise raise an exception
    1515"""
     16
     17from types import FloatType, IntType, LongType, NoneType
    1618
    1719from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
    … …  
    2325from anuga.caching import cache
    2426
    25 import Numeric as num
     27import anuga.utilities.numerical_tools as aunt
     28
     29import numpy as num
    2630
    2731
    … …  
    4347        if vertex_values is None:
    4448            N = len(domain)             # number_of_elements
    45             self.vertex_values = num.zeros((N, 3), num.Float)
     49            self.vertex_values = num.zeros((N, 3), num.float)
    4650        else:
    47             self.vertex_values = num.array(vertex_values, num.Float)
     51            self.vertex_values = num.array(vertex_values, num.float)
    4852
    4953            N, V = self.vertex_values.shape
    … …  
    5761
    5862        # Allocate space for other quantities
    59         self.centroid_values = num.zeros(N, num.Float)
    60         self.edge_values = num.zeros((N, 3), num.Float)
     63        self.centroid_values = num.zeros(N, num.float)
     64        self.edge_values = num.zeros((N, 3), num.float)
    6165
    6266        # Allocate space for Gradient
    63         self.x_gradient = num.zeros(N, num.Float)
    64         self.y_gradient = num.zeros(N, num.Float)
     67        self.x_gradient = num.zeros(N, num.float)
     68        self.y_gradient = num.zeros(N, num.float)
    6569
    6670        # Allocate space for Limiter Phi
    67         self.phi = num.zeros(N, num.Float)
     71        self.phi = num.zeros(N, num.float)
    6872
    6973        # Intialise centroid and edge_values
    … …  
    7276        # Allocate space for boundary values
    7377        L = len(domain.boundary)
    74         self.boundary_values = num.zeros(L, num.Float)
     78        self.boundary_values = num.zeros(L, num.float)
    7579
    7680        # Allocate space for updates of conserved quantities by
    … …  
    7882
    7983        # 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)
     84        self.explicit_update = num.zeros(N, num.float )
     85        self.semi_implicit_update = num.zeros(N, num.float )
     86        self.centroid_backup_values = num.zeros(N, num.float)
    8387
    8488        self.set_beta(1.0)
    … …  
    316320
    317321        numeric:
    318           Compatible list, Numeric array (see below) or constant.
     322          Compatible list, numeric array (see below) or constant.
    319323          If callable it will treated as a function (see below)
    320324          If instance of another Quantity it will be treated as such.
    … …  
    357361
    358362                  In case of location == 'centroids' the dimension values must
    359                   be a list of a Numerical array of length N,
     363                  be a list of a numerical array of length N,
    360364                  N being the number of elements.
    361365                  Otherwise it must be of dimension Nx3
    … …  
    402406
    403407        from anuga.geospatial_data.geospatial_data import Geospatial_data
    404         from types import FloatType, IntType, LongType, ListType, NoneType
    405408
    406409        # Treat special case: Polygon situation
    … …  
    457460            raise Exception, msg
    458461
    459         msg = 'Indices must be a list or None'
    460         assert type(indices) in [ListType, NoneType, num.ArrayType], msg
     462        msg = 'Indices must be a list, array or None'
     463        assert isinstance(indices, (NoneType, list, 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, (list, num.ndarray)):
    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             else:
    482                 msg = 'Illegal type for argument numeric: %s' % str(numeric)
    483                 raise msg
     481            else:   # see if it's coercible to a float (float, int or long, etc)
     482                try:
     483                    numeric = float(numeric)
     484                except ValueError:
     485                    msg = ("Illegal type for variable 'numeric': %s"
     486                           % type(numeric))
     487                    raise Exception, msg
     488                self.set_values_from_constant(numeric, location,
     489                                              indices, verbose)
    484490        elif quantity is not None:
    485491            self.set_values_from_quantity(quantity, location, indices, verbose)
    … …  
    515521            self.extrapolate_first_order()
    516522
    517 
    518 
    519523    ############################################################################
    520524    # Specific internal functions for setting values based on type
    … …  
    583587        """Set values for quantity
    584588
    585         values: Numeric array
     589        values: numeric array
    586590        location: Where values are to be stored.
    587591        Permissible options are: vertices, centroid, unique vertices
    … …  
    593597
    594598        In case of location == 'centroid' the dimension values must
    595         be a list of a Numerical array of length N, N being the number
     599        be a list of a numerical array of length N, N being the number
    596600        of elements.
    597601
    … …  
    607611        """
    608612
    609         values = num.array(values, num.Float)
     613        values = num.array(values, num.float)
    610614
    611615        if indices is not None:
    612             indices = num.array(indices, num.Int)
     616            indices = num.array(indices, num.int)
    613617            msg = ('Number of values must match number of indices: You '
    614618                   'specified %d values and %d indices'
    … …  
    637641                    'Values array must be 1d')
    638642
    639             self.set_vertex_values(values.flat, indices=indices,
     643            self.set_vertex_values(values.flatten(), indices=indices,
    640644                                   use_cache=use_cache, verbose=verbose)
    641645        else:
    … …  
    690694    # @param use_cache ??
    691695    # @param verbose True if this method is to be verbose.
    692     def set_values_from_function(self, f,
    693                                        location='vertices',
    694                                        indices=None,
    695                                        use_cache=False,
    696                                        verbose=False):
     696    def set_values_from_function(self,
     697                                 f,
     698                                 location='vertices',
     699                                 indices=None,
     700                                 use_cache=False,
     701                                 verbose=False):
    697702        """Set values for quantity using specified function
    698703
    … …  
    716721                indices = range(len(self))
    717722
    718             V = num.take(self.domain.get_centroid_coordinates(), indices)
     723            V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0)
    719724            x = V[:,0]; y = V[:,1]
    720725            if use_cache is True:
    … …  
    772777    # @param verbose ??
    773778    # @param use_cache ??
    774     def set_values_from_geospatial_data(self, geospatial_data,
    775                                               alpha,
    776                                               location,
    777                                               indices,
    778                                               verbose=False,
    779                                               use_cache=False):
     779    def set_values_from_geospatial_data(self,
     780                                        geospatial_data,
     781                                        alpha,
     782                                        location,
     783                                        indices,
     784                                        verbose=False,
     785                                        use_cache=False):
    780786        """Set values based on geo referenced geospatial data object."""
    781787
    … …  
    788794        from anuga.coordinate_transforms.geo_reference import Geo_reference
    789795
    790         points = ensure_numeric(points, num.Float)
    791         values = ensure_numeric(values, num.Float)
     796        points = ensure_numeric(points, num.float)
     797        values = ensure_numeric(values, num.float)
    792798
    793799        if location != 'vertices':
    … …  
    829835    # @param verbose True if this method is to be verbose.
    830836    # @param use_cache ??
    831     def set_values_from_points(self, points,
    832                                      values,
    833                                      alpha,
    834                                      location,
    835                                      indices,
    836                                      data_georef=None,
    837                                      verbose=False,
    838                                      use_cache=False):
     837    def set_values_from_points(self,
     838                               points,
     839                               values,
     840                               alpha,
     841                               location,
     842                               indices,
     843                               data_georef=None,
     844                               verbose=False,
     845                               use_cache=False):
    839846        """Set quantity values from arbitray data points using fit_interpolate.fit"""
    840847
    … …  
    851858    # @param use_cache
    852859    # @param max_read_lines
    853     def set_values_from_file(self, filename,
    854                                    attribute_name,
    855                                    alpha,
    856                                    location,
    857                                    indices,
    858                                    verbose=False,
    859                                    use_cache=False,
    860                                    max_read_lines=None):
     860    def set_values_from_file(self,
     861                             filename,
     862                             attribute_name,
     863                             alpha,
     864                             location,
     865                             indices,
     866                             verbose=False,
     867                             use_cache=False,
     868                             max_read_lines=None):
    861869        """Set quantity based on arbitrary points in a points file using
    862870        attribute_name selects name of attribute present in file.
    … …  
    10611069    # @param use_cache ??
    10621070    # @param verbose True if this method is to be verbose.
    1063     def get_interpolated_values(self, interpolation_points,
    1064                                       use_cache=False,
    1065                                       verbose=False):
     1071    def get_interpolated_values(self,
     1072                                interpolation_points,
     1073                                use_cache=False,
     1074                                verbose=False):
    10661075        """Get values at interpolation points
    10671076
    … …  
    10791088
    10801089        # Ensure that interpolation points is either a list of
    1081         # points, Nx2 array, or geospatial and convert to Numeric array
     1090        # points, Nx2 array, or geospatial and convert to numeric array
    10821091        if isinstance(interpolation_points, Geospatial_data):
    10831092            # Ensure interpolation points are in absolute UTM coordinates
    … …  
    11111120    # @param use_cache
    11121121    # @param verbose True if this method is to be verbose.
    1113     def get_values(self, 
     1122    def get_values(self,
    11141123                   interpolation_points=None,
    11151124                   location='vertices',
    … …  
    11191128        """Get values for quantity
    11201129
    1121         Extract values for quantity as a Numeric array.
     1130        Extract values for quantity as a numeric array.
    11221131
    11231132        Inputs:
    … …  
    11371146
    11381147        In case of location == 'centroids' the dimension of returned
    1139         values will be a list or a Numerical array of length N, N being
     1148        values will be a list or a numerical array of length N, N being
    11401149        the number of elements.
    11411150
    … …  
    11741183                            'edges', 'unique vertices']:
    11751184            msg = 'Invalid location: %s' % location
    1176             raise msg
     1185            raise Exception, msg
    11771186
    11781187        import types
    11791188
    1180         assert type(indices) in [types.ListType, types.NoneType,
    1181                                  num.ArrayType], \
    1182                    'Indices must be a list or None'
     1189        msg = "'indices' must be a list, array or None"
     1190        assert isinstance(indices, (NoneType, list, num.ndarray)), msg
    11831191
    11841192        if location == 'centroids':
    11851193            if (indices ==  None):
    11861194                indices = range(len(self))
    1187             return num.take(self.centroid_values, indices)
     1195            return num.take(self.centroid_values, indices, axis=0)
    11881196        elif location == 'edges':
    11891197            if (indices ==  None):
    11901198                indices = range(len(self))
    1191             return num.take(self.edge_values, indices)
     1199            return num.take(self.edge_values, indices, axis=0)
    11921200        elif location == 'unique vertices':
    11931201            if (indices ==  None):
    … …  
    12111219                    sum += self.vertex_values[triangle_id, vertex_id]
    12121220                vert_values.append(sum / len(triangles))
    1213             return num.array(vert_values, num.Float)
     1221            return num.array(vert_values, num.float)
    12141222        else:
    12151223            if (indices is None):
    12161224                indices = range(len(self))
    1217             return num.take(self.vertex_values, indices)
     1225            return num.take(self.vertex_values, indices, axis=0)
    12181226
    12191227    ##
    … …  
    12231231    # @param use_cache ??
    12241232    # @param verbose??
    1225     def set_vertex_values(self, A,
    1226                                 indices=None,
    1227                                 use_cache=False,
    1228                                 verbose=False):
     1233    def set_vertex_values(self,
     1234                          A,
     1235                          indices=None,
     1236                          use_cache=False,
     1237                          verbose=False):
    12291238        """Set vertex values for all unique vertices based on input array A
    12301239        which has one entry per unique vertex, i.e. one value for each row in
    … …  
    12371246
    12381247        # Check that A can be converted to array and is of appropriate dim
    1239         A = ensure_numeric(A, num.Float)
     1248        A = ensure_numeric(A, num.float)
    12401249        assert len(A.shape) == 1
    12411250
    … …  
    13371346
    13381347        if precision is None:
    1339             precision = num.Float
     1348            precision = num.float
    13401349
    13411350        if smooth is True:
    … …  
    13431352            V = self.domain.get_triangles()
    13441353            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    1345             A = num.zeros(N, num.Float)
     1354            A = num.zeros(N, num.float)
    13461355            points = self.domain.get_nodes()
    13471356
    … …  
    13611370                    if current_node == N:
    13621371                        msg = 'Current node exceeding number of nodes (%d) ' % N
    1363                         raise msg
     1372                        raise Exception, msg
    13641373
    13651374                    k += 1
    … …  
    13821391            V = self.domain.get_disconnected_triangles()
    13831392            points = self.domain.get_vertex_coordinates()
    1384             A = self.vertex_values.flat.astype(precision)
     1393            A = self.vertex_values.flatten().astype(precision)
    13851394
    13861395        # Return

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

    r7154 r7276  
    1111
    1212#include "Python.h"
    13 #include "Numeric/arrayobject.h"
     13#include "numpy/arrayobject.h"
    1414#include "math.h"
    1515
    … …  
    10281028        }
    10291029       
     1030        // check that numpy array objects arrays are C contiguous memory
     1031        CHECK_C_CONTIG(vertex_value_indices);
     1032        CHECK_C_CONTIG(number_of_triangles_per_node);
     1033        CHECK_C_CONTIG(vertex_values);
     1034        CHECK_C_CONTIG(A);
     1035
    10301036        N = vertex_value_indices -> dimensions[0];
    10311037        // printf("Got parameters, N=%d\n", N);

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

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

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

    r6717 r7276  
    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
    … …  
    549549                                      [0,0,9], [-6, 3, 3]])
    550550
    551         assert num.allclose(domain.quantities['stage'].centroid_values,
    552                             [2,5,3,0])
     551        assert num.allclose( domain.quantities['stage'].centroid_values,
     552                             [2,5,3,0] )
    553553
    554554        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
    … …  
    561561        domain.distribute_to_vertices_and_edges()
    562562
    563         # First order extrapolation
     563        #First order extrapolation
    564564        assert num.allclose( domain.quantities['stage'].vertex_values,
    565565                             [[ 2.,  2.,  2.],
    … …  
    614614
    615615        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    616         denom = num.ones(4, num.Float)-domain.timestep*sem
     616        denom = num.ones(4, num.float) - domain.timestep*sem
    617617
    618618#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
    … …  
    669669        domain.distribute_to_vertices_and_edges()
    670670
    671         # First order extrapolation
     671        #First order extrapolation
    672672        assert num.allclose( domain.quantities['stage'].vertex_values,
    673673                             [[ 2.,  2.,  2.],
    … …  
    766766       
    767767        #print domain.quantities['friction'].get_values()
     768        msg = ("domain.quantities['friction'].get_values()=\n%s\n"
     769               'should equal\n'
     770               '[[ 0.09,  0.09,  0.09],\n'
     771               ' [ 0.09,  0.09,  0.09],\n'
     772               ' [ 0.07,  0.07,  0.07],\n'
     773               ' [ 0.07,  0.07,  0.07],\n'
     774               ' [ 1.0,  1.0,  1.0],\n'
     775               ' [ 1.0,  1.0,  1.0]]'
     776               % str(domain.quantities['friction'].get_values()))
    768777        assert num.allclose(domain.quantities['friction'].get_values(),
    769778                            [[ 0.09,  0.09,  0.09],
    … …  
    772781                             [ 0.07,  0.07,  0.07],
    773782                             [ 1.0,  1.0,  1.0],
    774                              [ 1.0,  1.0,  1.0]])
     783                             [ 1.0,  1.0,  1.0]]), msg
    775784       
    776785        domain.set_region([set_bottom_friction, set_top_friction])
    … …  
    794803                             [ 11.0,  11.0,  11.0]])
    795804                             
    796 
    797805    def test_rectangular_periodic_and_ghosts(self):
    798806
    … …  
    909917
    910918#-------------------------------------------------------------
     919
    911920if __name__ == "__main__":
    912921    suite = unittest.makeSuite(Test_Domain,'test')

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

    r6145 r7276  
    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
    … …  
    183183       
    184184#-------------------------------------------------------------
     185
    185186if __name__ == "__main__":
    186187    suite = unittest.makeSuite(Test_ERMapper,'test')
    187188    runner = unittest.TextTestRunner()
    188189    runner.run(suite)
    189 
    190 
    191 
    192 

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

    r6717 r7276  
    1010from anuga.coordinate_transforms.geo_reference import Geo_reference
    1111
    12 import Numeric as num
     12import numpy as num
    1313
    1414
    … …  
    5353        e = [2.0, 2.0]
    5454        f = [4.0, 0.0]
    55 
    5655        nodes = num.array([a, b, c, d, e, f])
    5756
    5857        nodes_absolute = geo.get_absolute(nodes)
    5958       
    60         #bac, bce, ecf, dbe, daf, dae
    61         triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Int)
    62 
    63         domain = General_mesh(nodes, triangles,
    64                        geo_reference = geo)
    65         verts = domain.get_vertex_coordinates(triangle_id=0)       
    66         self.assert_(num.allclose(num.array([b,a,c]), verts))
     59        #                        bac,     bce,     ecf,     dbe
     60        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
     61
     62        domain = General_mesh(nodes, triangles, geo_reference=geo)
     63
     64        verts = domain.get_vertex_coordinates(triangle_id=0)    # bac
     65        msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s"
     66               % (str(num.array([b,a,c])), str(verts)))
     67        self.failUnless(num.allclose(num.array([b,a,c]), verts), msg)
     68
    6769        verts = domain.get_vertex_coordinates(triangle_id=0)       
    68         self.assert_(num.allclose(num.array([b,a,c]), verts))
     70        msg = ("num.array([b,a,c])=\n%s\nshould be close to 'verts'=\n%s"
     71               % (str(num.array([b,a,c])), str(verts)))
     72        self.assert_(num.allclose(num.array([b,a,c]), verts), msg)
     73
     74        verts = domain.get_vertex_coordinates(triangle_id=0, absolute=True)
     75        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
     76               % (str(num.array([nodes_absolute[1],
     77                                 nodes_absolute[0],
     78                                 nodes_absolute[2]])),
     79                  str(verts)))
     80        self.assert_(num.allclose(num.array([nodes_absolute[1],
     81                                             nodes_absolute[0],
     82                                             nodes_absolute[2]]),
     83                                  verts), msg)
     84
    6985        verts = domain.get_vertex_coordinates(triangle_id=0,
    7086                                              absolute=True)       
     87        msg = ("num.array([...])=\n%s\nshould be close to 'verts'=\n%s"
     88               % (str(num.array([nodes_absolute[1],
     89                                 nodes_absolute[0],
     90                                 nodes_absolute[2]])),
     91                  str(verts)))
    7192        self.assert_(num.allclose(num.array([nodes_absolute[1],
    72                                      nodes_absolute[0],
    73                                      nodes_absolute[2]]), verts))
    74         verts = domain.get_vertex_coordinates(triangle_id=0,
    75                                               absolute=True)       
    76         self.assert_(num.allclose(num.array([nodes_absolute[1],
    77                                      nodes_absolute[0],
    78                                      nodes_absolute[2]]), verts))
    79        
    80        
     93                                             nodes_absolute[0],
     94                                             nodes_absolute[2]]),
     95                                  verts), msg)
    8196
    8297    def test_get_vertex_coordinates_triangle_id(self):
    … …  
    117132        domain = General_mesh(nodes, triangles)
    118133
    119         value = [7]
    120         assert num.allclose(domain.get_triangles(), triangles)
     134        msg = ("domain.get_triangles()=\n%s\nshould be the same as "
     135               "'triangles'=\n%s"
     136               % (str(domain.get_triangles()), str(triangles)))
     137        assert num.allclose(domain.get_triangles(), triangles), msg
     138        msg = ("domain.get_triangles([0,4])=\n%s\nshould be the same as "
     139               "'[triangles[0], triangles[4]]' which is\n%s"
     140               % (str(domain.get_triangles([0,4])),
     141                  str([triangles[0], triangles[4]])))
    121142        assert num.allclose(domain.get_triangles([0,4]),
    122                         [triangles[0], triangles[4]])
     143                            [triangles[0], triangles[4]]), msg
    123144       
    124145
    … …  
    271292        nodes_absolute = geo.get_absolute(nodes)
    272293       
    273         #bac, bce, ecf, dbe, daf, dae
     294        #                        bac,     bce,     ecf,     dbe
    274295        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    275296
    276         domain = General_mesh(nodes, triangles,
    277                        geo_reference = geo)
     297        domain = General_mesh(nodes, triangles, geo_reference = geo)
    278298        node = domain.get_node(2)       
    279         self.assertEqual(c, node)
     299        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
     300        self.failUnless(num.alltrue(c == node), msg)
     301
     302        # repeat get_node(), see if result same
     303        node = domain.get_node(2)       
     304        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
     305        self.failUnless(num.alltrue(c == node), msg)
    280306       
    281307        node = domain.get_node(2, absolute=True)     
    282         self.assertEqual(nodes_absolute[2], node)
    283        
     308        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
     309               % (str(nodes_absolute[2]), str(node)))
     310        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
     311       
     312        # repeat get_node(2, absolute=True), see if result same
    284313        node = domain.get_node(2, absolute=True)     
    285         self.assertEqual(nodes_absolute[2], node)
     314        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
     315               % (str(nodes_absolute[2]), str(node)))
     316        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
    286317       
    287318
    … …  
    321352        self.failUnlessRaises(AssertionError, General_mesh,
    322353                              nodes, triangles, geo_reference=geo)
    323        
    324  
    325 
    326 #-------------------------------------------------------------
     354
     355################################################################################
     356
    327357if __name__ == "__main__":
    328     suite = unittest.makeSuite(Test_General_Mesh,'test')
    329     #suite = unittest.makeSuite(Test_General_Mesh,'test_get_node')   
     358    suite = unittest.makeSuite(Test_General_Mesh, 'test')
    330359    runner = unittest.TextTestRunner()
    331360    runner.run(suite)

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

    r6195 r7276  
    88from anuga.config import epsilon
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    270270
    271271
    272 
    273 
    274 
    275272#-------------------------------------------------------------
     273
    276274if __name__ == "__main__":
    277275    suite = unittest.makeSuite(Test_Generic_Boundary_Conditions, 'test')
    278276    runner = unittest.TextTestRunner()
    279277    runner.run(suite)
    280 
    281 
    282 

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

    r6195 r7276  
    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.
    43 
    44 
     44        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    4545
    4646
    4747#-------------------------------------------------------------
     48
    4849if __name__ == "__main__":
    4950    suite = unittest.makeSuite(Test_Domain,'test')

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

    r6717 r7276  
    1919from anuga.utilities.numerical_tools import ensure_numeric
    2020
    21 import Numeric as num
     21import numpy as num
    2222
    2323
    … …  
    993993                  [  75735.4765625 ,  23762.00585938],
    994994                  [  52341.70703125,  38563.39453125]]
    995 
    996         ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
    997995
    998996        triangles = [[19, 0,15],
    … …  
    10971095                  [  35406.3359375 ,  79332.9140625 ]]
    10981096
    1099         scaled_points = ensure_numeric(points, num.Int)/1000  # Simplify for ease of interpretation
     1097        scaled_points = ensure_numeric(points, num.int)/1000  # Simplify for ease of interpretation
    11001098
    11011099        triangles = [[ 0, 1, 2],
    … …  
    11761174        """
    11771175       
    1178        
    11791176        # test
    11801177        a = [0.0, 0.0]
    … …  
    11831180
    11841181        absolute_points = [a, b, c]
    1185         vertices = [[0,1,2]]
    1186        
    1187         geo = Geo_reference(56,67,-56)
     1182        vertices = [[0, 1, 2]]
     1183       
     1184        geo = Geo_reference(56, 67, -56)
    11881185
    11891186        relative_points = geo.change_points_geo_ref(absolute_points)
    1190 
    1191         #print 'Relative', relative_points
    1192         #print 'Absolute', absolute_points       
    11931187
    11941188        mesh = Mesh(relative_points, vertices, geo_reference=geo)
    … …  
    18461840
    18471841
    1848 
    1849        
    18501842#-------------------------------------------------------------
     1843
    18511844if __name__ == "__main__":
    1852     #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
    1853     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_partially_coinciding')
    1854     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments7')
    1855     suite = unittest.makeSuite(Test_Mesh,'test')
     1845    suite = unittest.makeSuite(Test_Mesh, 'test')
    18561846    runner = unittest.TextTestRunner()
    18571847    runner.run(suite)
    1858 
    1859 
    1860 
    1861 

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

    r6191 r7276  
    1414from anuga.pmesh.mesh import importMeshFromFile
    1515
    16 import Numeric as num
     16import numpy as num
    1717
    1818
    … …  
    255255
    256256#-------------------------------------------------------------
     257
    257258if __name__ == "__main__":
    258259    suite = unittest.makeSuite(Test_pmesh2domain, 'test')
    259260    runner = unittest.TextTestRunner()
    260261    runner.run(suite)
    261 
    262 
    263 
    264 

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

    r6195 r7276  
    1515from anuga.utilities.polygon import *
    1616
    17 import Numeric as num
     17import numpy as num
    1818
    1919
    … …  
    504504        quantity.set_values(0.0)
    505505        quantity.set_values(3.14, polygon=polygon)
     506        msg = ('quantity.vertex_values=\n%s\nshould be close to\n'
     507               '[[0,0,0],\n'
     508               ' [3.14,3.14,3.14],\n'
     509               ' [3.14,3.14,3.14],\n'
     510               ' [0,0,0]]' % str(quantity.vertex_values))
    506511        assert num.allclose(quantity.vertex_values,
    507512                            [[0,0,0],
    508513                             [3.14,3.14,3.14],
    509514                             [3.14,3.14,3.14],                         
    510                              [0,0,0]])               
     515                             [0,0,0]]), msg
    511516
    512517
    … …  
    17651770        quantity = Quantity(self.mesh4)
    17661771
    1767         quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Float)
     1772        quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.float)
    17681773
    17691774        quantity.interpolate_from_vertices_to_edges()
    … …  
    17811786                                          [3., 2.5, 1.5],
    17821787                                          [3.5, 4.5, 3.],
    1783                                           [2.5, 3.5, 2]], num.Float)
     1788                                          [2.5, 3.5, 2]], num.float)
    17841789
    17851790        quantity.interpolate_from_edges_to_vertices()
    … …  
    18351840
    18361841        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1837         denom = num.ones(4, num.Float)-timestep*sem
     1842        denom = num.ones(4, num.float)-timestep*sem
    18381843
    18391844        x = num.array([1, 2, 3, 4])/denom
    … …  
    18591864
    18601865        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1861         denom = num.ones(4, num.Float)-timestep*sem
     1866        denom = num.ones(4, num.float)-timestep*sem
    18621867
    18631868        x = num.array([1., 2., 3., 4.])
    … …  
    19011906
    19021907        bed = domain.quantities['elevation'].vertex_values
    1903         stage = num.zeros(bed.shape, num.Float)
     1908        stage = num.zeros(bed.shape, num.float)
    19041909
    19051910        h = 0.03
    … …  
    19851990
    19861991        bed = domain.quantities['elevation'].vertex_values
    1987         stage = num.zeros(bed.shape, num.Float)
     1992        stage = num.zeros(bed.shape, num.float)
    19881993
    19891994        h = 0.03
    … …  
    19992004        stage = domain.quantities['stage']
    20002005        A, V = stage.get_vertex_values(xy=False, smooth=False)
    2001         Q = stage.vertex_values.flat
     2006        Q = stage.vertex_values.flatten()
    20022007
    20032008        for k in range(8):
    … …  
    25082513
    25092514#-------------------------------------------------------------
     2515
    25102516if __name__ == "__main__":
    25112517    suite = unittest.makeSuite(Test_Quantity, 'test')   
    2512     #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')
    25172518    runner = unittest.TextTestRunner()
    25182519    runner.run(suite)

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_region.py

    r6145 r7276  
    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],
    … …  
    189189                         
    190190    def test_unique_vertices_average_loc_vert(self):
    191         """
    192         get values based on triangle lists.
    193         """
    194         from mesh_factory import rectangular
    195         from shallow_water import Domain
    196 
    197         #Create basic mesh
    198         points, vertices, boundary = rectangular(1, 3)
    199         #Create shallow water domain
    200         domain = Domain(points, vertices, boundary)
    201         domain.build_tagged_elements_dictionary({'bottom':[0,1],
    202                                                  'top':[4,5],
    203                                                  'not_bottom':[2,3,4,5]})
     191        """Get values based on triangle lists."""
     192
     193        from mesh_factory import rectangular
     194        from shallow_water import Domain
     195
     196        #Create basic mesh
     197        points, vertices, boundary = rectangular(1, 3)
     198
     199        #Create shallow water domain
     200        domain = Domain(points, vertices, boundary)
     201        domain.build_tagged_elements_dictionary({'bottom': [0, 1],
     202                                                 'top': [4, 5],
     203                                                 'not_bottom': [2, 3, 4, 5]})
    204204
    205205        #Set friction
    206206        domain.set_quantity('friction', add_x_y)
    207         av_bottom = 2.0/3.0
     207        av_bottom = 2.0 / 3.0
    208208        add = 60.0
    209209        calc_frict = av_bottom + add
    210210        domain.set_region(Add_value_to_region('bottom', 'friction', add,
    211                           initial_quantity='friction',
    212                            #location='unique vertices',
    213                            location='vertices',
    214                            average=True
    215                           ))
    216 
    217         #print domain.quantities['friction'].get_values()
     211                                              initial_quantity='friction',
     212                                              location='vertices',
     213                                              average=True))
     214
    218215        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])
     216
     217        expected = [calc_frict, calc_frict, calc_frict]
     218        msg = ('frict_points[0]=%s\nexpected=%s'
     219               % (str(frict_points[0]), str(expected)))
     220        assert num.allclose(frict_points[0], expected), msg
     221
     222        msg = ('frict_points[1]=%s\nexpected=%s'
     223               % (str(frict_points[1]), str(expected)))
     224        assert num.allclose(frict_points[1], expected), msg
    223225 
    224226    def test_unique_vertices_average_loc_unique_vert(self):
    … …  
    261263                         
    262264#-------------------------------------------------------------
     265
    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)

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r6318 r7276  
    2222import string
    2323
    24 import Numeric as num
     24import numpy as num
    2525
    2626
    … …  
    183183
    184184        last_time_index = len(time)-1 #Last last_time_index
    185         d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1))
    186         d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
    187         d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
     185        d_stage = num.reshape(num.take(stage[last_time_index, :],
     186                                       [0,5,10,15],
     187                                       axis=0),
     188                              (4,1))
     189        d_uh = num.reshape(num.take(xmomentum[last_time_index, :],
     190                                    [0,5,10,15],
     191                                   axis=0),
     192                           (4,1))
     193        d_vh = num.reshape(num.take(ymomentum[last_time_index, :],
     194                                    [0,5,10,15],
     195                                   axis=0),
     196                           (4,1))
    188197        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    189198
    … …  
    195204
    196205        #And the midpoints are found now
    197         Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15])
    198         Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15])
     206        Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0)
     207        Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0)
    199208
    200209        diag = num.concatenate( (Dx, Dy), axis=1)
    … …  
    219228
    220229        timestep = 0 #First timestep
    221         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    222         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    223         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     230        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     231        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     232        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    224233        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    225234
    … …  
    241250
    242251        timestep = 33
    243         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    244         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    245         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     252        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     253        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     254        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    246255        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    247256
    … …  
    262271
    263272        timestep = 15
    264         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    265         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    266         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     273        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     274        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     275        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    267276        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    268277
    … …  
    275284        #
    276285        timestep = 16
    277         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    278         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    279         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     286        d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15], axis=0), (4,1))
     287        d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
     288        d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15], axis=0), (4,1))
    280289        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    281290
    … …  
    377386        x = fid.variables['x'][:]
    378387        y = fid.variables['y'][:]
     388        # we 'cast' to 64 bit floats to pass this test
     389        # SWW file quantities are stored as 32 bits
     390        x = num.array(x, num.float)
     391        y = num.array(y, num.float)
     392
    379393        stage = fid.variables['stage'][:]
    380394        xmomentum = fid.variables['xmomentum'][:]
    … …  
    386400
    387401        last_time_index = len(time)-1 #Last last_time_index     
    388         d_stage = num.reshape(num.take(stage[last_time_index, :], [0,5,10,15]), (4,1))
    389         d_uh = num.reshape(num.take(xmomentum[last_time_index, :], [0,5,10,15]), (4,1))
    390         d_vh = num.reshape(num.take(ymomentum[last_time_index, :], [0,5,10,15]), (4,1))
     402        d_stage = num.reshape(num.take(stage[last_time_index, :],
     403                                       [0,5,10,15],
     404                                       axis=0),
     405                              (4,1))
     406        d_uh = num.reshape(num.take(xmomentum[last_time_index, :],
     407                                    [0,5,10,15],
     408                                   axis=0),
     409                           (4,1))
     410        d_vh = num.reshape(num.take(ymomentum[last_time_index, :],
     411                                    [0,5,10,15],
     412                                    axis=0),
     413                           (4,1))
    391414        D = num.concatenate((d_stage, d_uh, d_vh), axis=1)
    392415
    … …  
    398421
    399422        #And the midpoints are found now
    400         Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15])
    401         Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15])
    402 
    403         diag = num.concatenate( (Dx, Dy), axis=1)
     423        Dx = num.take(num.reshape(x, (16,1)), [0,5,10,15], axis=0)
     424        Dy = num.take(num.reshape(y, (16,1)), [0,5,10,15], axis=0)
     425
     426        diag = num.concatenate((Dx, Dy), axis=1)
    404427        d_midpoints = (diag[1:] + diag[:-1])/2
    405428
    … …  
    415438
    416439        t = time[last_time_index]                         
    417         q = f(t, point_id=0); assert num.allclose(r0, q)
    418         q = f(t, point_id=1); assert num.allclose(r1, q)
    419         q = f(t, point_id=2); assert num.allclose(r2, q)
     440
     441        q = f(t, point_id=0)
     442        msg = '\nr0=%s\nq=%s' % (str(r0), str(q))
     443        assert num.allclose(r0, q), msg
     444
     445        q = f(t, point_id=1)
     446        msg = '\nr1=%s\nq=%s' % (str(r1), str(q))
     447        assert num.allclose(r1, q), msg
     448
     449        q = f(t, point_id=2)
     450        msg = '\nr2=%s\nq=%s' % (str(r2), str(q))
     451        assert num.allclose(r2, q), msg
    420452
    421453
    … …  
    424456
    425457        timestep = 0 #First timestep
    426         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    427         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    428         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     458        d_stage = num.reshape(num.take(stage[timestep, :],
     459                                       [0,5,10,15],
     460                                       axis=0),
     461                              (4,1))
     462        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     463                                    [0,5,10,15],
     464                                    axis=0),
     465                           (4,1))
     466        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     467                                    [0,5,10,15],
     468                                    axis=0),
     469                           (4,1))
    429470        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    430471
    … …  
    446487
    447488        timestep = 33
    448         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    449         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    450         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     489        d_stage = num.reshape(num.take(stage[timestep, :],
     490                                       [0,5,10,15],
     491                                       axis=0),
     492                              (4,1))
     493        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     494                                    [0,5,10,15],
     495                                    axis=0),
     496                           (4,1))
     497        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     498                                    [0,5,10,15],
     499                                    axis=0),
     500                           (4,1))
    451501        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    452502
    … …  
    467517
    468518        timestep = 15
    469         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    470         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    471         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     519        d_stage = num.reshape(num.take(stage[timestep, :],
     520                                       [0,5,10,15],
     521                                       axis=0),
     522                              (4,1))
     523        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     524                                    [0,5,10,15],
     525                                    axis=0),
     526                           (4,1))
     527        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     528                                    [0,5,10,15],
     529                                    axis=0),
     530                           (4,1))
    472531        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    473532
    … …  
    480539        #
    481540        timestep = 16
    482         d_stage = num.reshape(num.take(stage[timestep, :], [0,5,10,15]), (4,1))
    483         d_uh = num.reshape(num.take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
    484         d_vh = num.reshape(num.take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
     541        d_stage = num.reshape(num.take(stage[timestep, :],
     542                                       [0,5,10,15],
     543                                       axis=0),
     544                              (4,1))
     545        d_uh = num.reshape(num.take(xmomentum[timestep, :],
     546                                    [0,5,10,15],
     547                                    axis=0),
     548                           (4,1))
     549        d_vh = num.reshape(num.take(ymomentum[timestep, :],
     550                                    [0,5,10,15],
     551                                    axis=0),
     552                           (4,1))
    485553        D = num.concatenate( (d_stage, d_uh, d_vh), axis=1)
    486554
    … …  
    630698                    q1 = F(t+60, point_id=id)
    631699
    632                 if q0 == NAN:
     700                if num.alltrue(q0 == NAN):
    633701                    actual = q0
    634702                else:
    … …  
    641709                #print "actual", actual
    642710                #print
    643                 if q0 == NAN:
    644                      self.failUnless( q == actual, 'Fail!')
     711                if num.alltrue(q0 == NAN):
     712                     self.failUnless(num.alltrue(q == actual), 'Fail!')
    645713                else:
    646714                    assert num.allclose(q, actual)
    … …  
    11811249        #FIXME: Division is not expected to work for integers.
    11821250        #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)                 
     1251        foo = num.array([[1,2,3], [4,5,6]], num.float)
     1252
     1253        bar = num.array([[-1,0,5], [6,1,1]], num.float)                 
    11861254
    11871255        D = {'X': foo, 'Y': bar}
    … …  
    12011269
    12021270        # make an error for zero on zero
    1203         # this is really an error in Numeric, SciPy core can handle it
     1271        # this is really an error in numeric, SciPy core can handle it
    12041272        # Z = apply_expression_to_dictionary('0/Y', D)
    12051273
    … …  
    19191987        if 314 < angle < 316: v=1
    19201988        assert v==1
    1921 
    1922  
    1923 
    1924 
    19251989       
    19261990
    19271991#-------------------------------------------------------------
     1992
    19281993if __name__ == "__main__":
    1929     suite = unittest.makeSuite(Test_Util,'test')
    1930 #    suite = unittest.makeSuite(Test_Util,'test_sww2csv_gauges')
     1994    suite = unittest.makeSuite(Test_Util, 'test')
    19311995#    runner = unittest.TextTestRunner(verbosity=2)
    19321996    runner = unittest.TextTestRunner(verbosity=1)
    19331997    runner.run(suite)
    1934 
    1935 
    1936 
    1937 

  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r7012 r7276  
    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
    … …  
    476476        if boundary_polygon is not None:
    477477            #removes sts points that do not lie on boundary
    478             quantities[name] = num.take(quantities[name], gauge_id, 1)
     478            quantities[name] = num.take(quantities[name], gauge_id, axis=1)
    479479           
    480480    # Close sww, tms or sts netcdf file         
    … …  
    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                     
    … …  
    12561256    n0 = int(n0)
    12571257    m = len(locations)
    1258     model_time = num.zeros((n0, m, p), num.Float)
    1259     stages = num.zeros((n0, m, p), num.Float)
    1260     elevations = num.zeros((n0, m, p), num.Float)
    1261     momenta = num.zeros((n0, m, p), num.Float)
    1262     xmom = num.zeros((n0, m, p), num.Float)
    1263     ymom = num.zeros((n0, m, p), num.Float)
    1264     speed = num.zeros((n0, m, p), num.Float)
    1265     bearings = num.zeros((n0, m, p), num.Float)
    1266     due_east = 90.0*num.ones((n0, 1), num.Float)
    1267     due_west = 270.0*num.ones((n0, 1), num.Float)
    1268     depths = num.zeros((n0, m, p), num.Float)
    1269     eastings = num.zeros((n0, m, p), num.Float)
     1258    model_time = num.zeros((n0, m, p), num.float)
     1259    stages = num.zeros((n0, m, p), num.float)
     1260    elevations = num.zeros((n0, m, p), num.float)
     1261    momenta = num.zeros((n0, m, p), num.float)
     1262    xmom = num.zeros((n0, m, p), num.float)
     1263    ymom = num.zeros((n0, m, p), num.float)
     1264    speed = num.zeros((n0, m, p), num.float)
     1265    bearings = num.zeros((n0, m, p), num.float)
     1266    due_east = 90.0*num.ones((n0, 1), num.float)
     1267    due_west = 270.0*num.ones((n0, 1), num.float)
     1268    depths = num.zeros((n0, m, p), num.float)
     1269    eastings = num.zeros((n0, m, p), num.float)
    12701270    min_stages = []
    12711271    max_stages = []
    … …  
    12791279    min_speeds = []   
    12801280    max_depths = []
    1281     model_time_plot3d = num.zeros((n0, m), num.Float)
    1282     stages_plot3d = num.zeros((n0, m), num.Float)
    1283     eastings_plot3d = num.zeros((n0, m),num.Float)
     1281    model_time_plot3d = num.zeros((n0, m), num.float)
     1282    stages_plot3d = num.zeros((n0, m), num.float)
     1283    eastings_plot3d = num.zeros((n0, m),num.float)
    12841284    if time_unit is 'mins': scale = 60.0
    12851285    if time_unit is 'hours': scale = 3600.0
    … …  
    18001800        # Remove the loners from verts
    18011801        # Could've used X=compress(less(loners,N),loners)
    1802         # verts=num.take(verts,X)  to Remove the loners from verts
     1802        # verts=num.take(verts,X,axis=0)  to Remove the loners from verts
    18031803        # but I think it would use more memory
    18041804        new_i = lone_start      # point at first loner - 'shuffle down' target
    … …  
    18341834    """
    18351835       
    1836     xc = num.zeros(triangles.shape[0], num.Float) # Space for centroid info
     1836    xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info
    18371837   
    18381838    for k in range(triangles.shape[0]):
    … …  
    21212121                #add tide to stage if provided
    21222122                if quantity == 'stage':
    2123                      quantity_value[quantity] = num.array(quantity_value[quantity], num.Float) \
    2124                                                           + directory_add_tide
     2123                     quantity_value[quantity] = num.array(quantity_value[quantity],
     2124                                                          num.float) + directory_add_tide
    21252125
    21262126                #condition to find max and mins for all the plots
    … …  
    23702370        file.close()
    23712371
    2372 
    23732372##
    23742373# @brief ??
    … …  
    23892388   
    23902389    Inputs:
    2391        
    23922390        NOTE: if using csv2timeseries_graphs after creating csv file,
    23932391        it is essential to export quantities 'depth' and 'elevation'.
    … …  
    24142412           myfile_2_point1.csv if <out_name> ='myfile_2_'
    24152413           
    2416            
    24172414        They will all have a header
    24182415   
    … …  
    24382435    import string
    24392436    from anuga.shallow_water.data_manager import get_all_swwfiles
    2440 
    2441 #    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
    2442     #print "points",points
    24432437
    24442438    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
    … …  
    24572451    point_name = []
    24582452   
    2459     #read point info from file
     2453    # read point info from file
    24602454    for i,row in enumerate(point_reader):
    2461         #read header and determine the column numbers to read correcty.
     2455        # read header and determine the column numbers to read correctly.
    24622456        if i==0:
    24632457            for j,value in enumerate(row):
    … …  
    24712465       
    24722466    #convert to array for file_function
    2473     points_array = num.array(points,num.Float)
     2467    points_array = num.array(points,num.float)
    24742468       
    24752469    points_array = ensure_absolute(points_array)
    … …  
    25252519    for sww_file in sww_files:
    25262520        sww_file = join(dir_name, sww_file+'.sww')
    2527         #print 'sww file = ',sww_file
    25282521        callable_sww = file_function(sww_file,
    25292522                                     quantities=core_quantities,
    … …  
    25792572                                    momentum = sqrt(point_quantities[2]**2
    25802573                                                    + point_quantities[3]**2)
    2581         #                            vel = momentum/depth             
    25822574                                    vel = momentum / (point_quantities[0]
    25832575                                                      - point_quantities[1])
    2584         #                            vel = momentum/(depth + 1.e-6/depth)
    25852576                                else:
    25862577                                    momentum = 0
    … …  
    25932584                                                            point_quantities[3]))
    25942585
    2595                 #print 'point list before write (writer %s) = %s' % (str(point_name[point_i]), str(points_list))
    25962586                points_writer[point_i].writerow(points_list)
    2597            
    25982587
    25992588##

  • anuga_core/source/anuga/advection/advection.py

    r6146 r7276  
    3333from anuga.abstract_2d_finite_volumes.domain import *
    3434
    35 import Numeric as num
     35import numpy as num
    3636
    3737
    … …  
    252252        stage_bdry = Stage.boundary_values
    253253
    254         flux = num.zeros(1, num.Float) #Work array for summing up fluxes
     254        flux = num.zeros(1, num.float) #Work array for summing up fluxes
    255255
    256256        #Loop

  • anuga_core/source/anuga/advection/advection_ext.c

    r5162 r7276  
    1010
    1111#include "Python.h"
    12 #include "Numeric/arrayobject.h"
     12#include "numpy/arrayobject.h"
    1313#include "math.h"
    1414#include "stdio.h"

  • anuga_core/source/anuga/advection/test_advection.py

    r6146 r7276  
    88from anuga.advection.advection import Domain, Transmissive_boundary, Dirichlet_boundary
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212

  • anuga_core/source/anuga/alpha_shape/alpha_shape.py

    r6268 r7276  
    2525from anuga.geospatial_data.geospatial_data import Geospatial_data
    2626
    27 import Numeric as num
     27import numpy as num
    2828
    2929
    … …  
    8888                raise PointError, "Three points on a straight line"
    8989       
    90         #Convert input to Numeric arrays
    91         self.points = num.array(points, num.Float)
     90        #Convert input to numeric arrays
     91        self.points = num.array(points, num.float)
    9292
    9393   
    … …  
    289289                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
    290290
    291         if num.alltrue(denom != 0.0):               
    292             dx = num.divide_safe(y31*dist21 - y21*dist31,denom)
    293             dy = num.divide_safe(x21*dist31 - x31*dist21,denom)
    294         else:
    295             raise  AlphaError
    296            
     291        if num.any(denom == 0.0):
     292            raise AlphaError
     293
     294        dx = num.divide(y31*dist21 - y21*dist31, denom)
     295        dy = num.divide(x21*dist31 - x31*dist21, denom)
     296
    297297        self.triradius = 0.5*num.sqrt(dx*dx + dy*dy)
    298298        #print "triangle radii", self.triradius

  • anuga_core/source/anuga/alpha_shape/test_alpha_shape.py

    r6147 r7276  
    55import unittest
    66
    7 import Numeric as num
     7import numpy as num
    88
    99try:
    … …  
    394394                  (46, 45)]
    395395        assert num.allclose(answer, result)   
     396
    396397#-------------------------------------------------------------
     398
    397399if __name__ == "__main__":
    398     #print "starting tests \n"
    399400    suite = unittest.makeSuite(TestCase,'test')
    400401    runner = unittest.TextTestRunner(verbosity=1)
    401402    runner.run(suite)
    402     #print "finished tests \n"
    403    
    404 
    405    
    406    
    407 
    408 
    409 

  • anuga_core/source/anuga/caching/caching.py

    r7197 r7276  
    5151  unix = True
    5252
    53 import Numeric as num
     53import numpy as num
    5454
    5555
    … …  
    13941394      I.sort()   
    13951395      val = myhash(I, ids)
    1396   elif type(T) == num.ArrayType:
     1396  elif isinstance(T, num.ndarray):
    13971397      T = num.array(T) # Ensure array is contiguous
    13981398
    … …  
    14651465            identical = compare(a, b, ids)
    14661466           
    1467     elif type(A) == num.ArrayType:
     1467    elif isinstance(A, num.ndarray):
    14681468        # Use element by element comparison
    14691469        identical = num.alltrue(A==B)
    … …  
    24392439      argstr = argstr + "'"+str(args)+"'"
    24402440    else:
    2441       # Truncate large Numeric arrays before using str()
    2442       if type(args) == num.ArrayType:
     2441      # Truncate large numeric arrays before using str()
     2442      if isinstance(args, num.ndarray):
    24432443#        if len(args.flat) > textwidth: 
    24442444#        Changed by Duncan and Nick 21/2/07 .flat has problems with

  • anuga_core/source/anuga/caching/test_caching.py

    r6406 r7276  
    77from anuga.caching.dummy_classes_for_testing import Dummy, Dummy_memorytest
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    2626 
    2727def f_numeric(A, B):
    28   """Operation on Numeric arrays
     28  """Operation on numeric arrays
    2929  """
    3030 
    … …  
    123123        """test_caching_of_numeric_arrays
    124124       
    125         Test that Numeric arrays can be recognised by caching even if their id's are different
     125        Test that numeric arrays can be recognised by caching even if their id's are different
    126126        """
    127127       
    … …  
    160160
    161161
    162             assert T1 == T2, 'Cached result does not match computed result'
    163             assert T2 == T3, 'Cached result does not match computed result'
     162            assert num.alltrue(T1 == T2), 'Cached result does not match computed result'
     163            assert num.alltrue(T2 == T3), 'Cached result does not match computed result'
    164164           
    165165
    166166    def test_hash_collision(self):
    167         """test_hash_collision(self):
    168        
    169         Test that hash collisons are dealt with correctly
    170         """
     167        """Test that hash collisons are dealt with correctly"""
    171168       
    172169        verbose = False
    … …  
    202199            T2 = cache(f_numeric, (A1, A1),
    203200                       compression=comp, verbose=verbose)
    204            
     201           
     202            T1_ref = f_numeric(A0, A0)
     203            T2_ref = f_numeric(A1, A1)
    205204            T1_ref = f_numeric(A0, A0)
    206205            T2_ref = f_numeric(A1, A1)
    … …  
    419418
    420419       
    421         x = num.arange(10).astype(num.Float)
     420        x = num.arange(10).astype(num.float)
    422421       
    423422        ref1 = f1(x)
    … …  
    903902#-------------------------------------------------------------
    904903if __name__ == "__main__":
    905     #suite = unittest.makeSuite(Test_Caching, 'test_caching_of_callable_objects')
    906904    suite = unittest.makeSuite(Test_Caching, 'test')
    907905    runner = unittest.TextTestRunner()

  • anuga_core/source/anuga/config.py

    r7055 r7276  
    44import os
    55import sys
     6
    67
    78################################################################################
    … …  
    162163                             # Too large (100) creates 'flopping' water
    163164                             # Too small (0) creates 'creep'
    164                            
     165
    165166maximum_froude_number = 100.0 # To be used in limiters.
    166167
    … …  
    183184
    184185################################################################################
     186# NetCDF-specific type constants.  Used when defining NetCDF file variables.
     187################################################################################
     188
     189netcdf_char = 'c'
     190netcdf_byte = 'b'
     191netcdf_int = 'i'
     192netcdf_float = 'd'
     193netcdf_float64 = 'd'
     194netcdf_float32 = 'f'
     195
     196################################################################################