Changeset 6428


Ignore:
Timestamp:
Feb 27, 2009, 11:54:09 AM (16 years ago)
Author:
rwilson
Message:

numpy changes

Location:
branches/numpy/anuga
Files:
16 edited

Legend:

Unmodified
Added
Removed

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

    r6410 r6428  
    3737        triangles = [ [1,0,2], [1,2,3] ]   # bac, bce
    3838
    39         # Create mesh with two triangles: bac and bce   
     39        # Create mesh with two triangles: bac and bce
    4040        mesh = Mesh(nodes, triangles)
    4141
    … …  
    5656    """
    5757
    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,                 
     58    # FIXME: It would be a good idea to use geospatial data as an alternative
     59    #        input
     60    def __init__(self,
     61                 nodes,
     62                 triangles,
     63                 geo_reference=None,
    6264                 number_of_full_nodes=None,
    63                  number_of_full_triangles=None,                 
     65                 number_of_full_triangles=None,
    6466                 verbose=False):
    6567        """Build triangular 2d mesh from nodes and triangle information
    6668
    6769        Input:
    68        
     70
    6971          nodes: x,y coordinates represented as a sequence of 2-tuples or
    7072                 a Nx2 numeric array of floats.
    71                  
     73
    7274          triangles: sequence of 3-tuples or Mx3 numeric array of
    7375                     non-negative integers representing indices into
    7476                     the nodes array.
    75        
     77
    7678          georeference (optional): If specified coordinates are
    7779          assumed to be relative to this origin.
    … …  
    8385        In this case it is usefull to specify the number of real (called full)
    8486        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' 
     87
     88        """
     89
     90        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
    8991
    9092        self.triangles = num.array(triangles, num.int)
    9193        self.nodes = num.array(nodes, num.float)
    9294
    93 
    94         # Register number of elements and nodes
     95        # Register number of elements and nodes
    9596        self.number_of_triangles = N = self.triangles.shape[0]
    96         self.number_of_nodes = self.nodes.shape[0]       
    97 
    98        
     97        self.number_of_nodes = self.nodes.shape[0]
    9998
    10099        if number_of_full_nodes is None:
    … …  
    102101        else:
    103102            assert int(number_of_full_nodes)
    104             self.number_of_full_nodes = number_of_full_nodes           
    105 
     103            self.number_of_full_nodes = number_of_full_nodes
    106104
    107105        if number_of_full_triangles is None:
    108             self.number_of_full_triangles = self.number_of_triangles           
     106            self.number_of_full_triangles = self.number_of_triangles
    109107        else:
    110             assert int(number_of_full_triangles)           
     108            assert int(number_of_full_triangles)
    111109            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            
    118110
    119111        # FIXME: this stores a geo_reference, but when coords are returned
    120112        # This geo_ref is not taken into account!
    121113        if geo_reference is None:
    122             self.geo_reference = Geo_reference() #Use defaults
     114            self.geo_reference = Geo_reference()    # Use defaults
    123115        else:
    124116            self.geo_reference = geo_reference
    125117
    126118        # 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)
     119        msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
     120               'The supplied array has the shape: %s'
     121               % str(self.triangles.shape))
    130122        assert len(self.triangles.shape) == 2, msg
    131123
    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)
     124        msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples'
     125               'The supplied array has the shape: %s' % str(self.nodes.shape))
    135126        assert len(self.nodes.shape) == 2, msg
    136127
    … …  
    138129        assert max(self.triangles.flat) < self.nodes.shape[0], msg
    139130
    140 
    141131        # FIXME: Maybe move to statistics?
    142132        # 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]) ]
     133        xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]),
     134                     max(self.nodes[:,0]), max(self.nodes[:,1])]
    145135
    146136        self.xy_extent = num.array(xy_extent, num.float)
    147 
    148137
    149138        # Allocate space for geometric quantities
    … …  
    155144        self.vertex_coordinates = V = self.compute_vertex_coordinates()
    156145
    157 
    158146        # Initialise each triangle
    159147        if verbose:
    160             print 'General_mesh: Computing areas, normals and edgelenghts'
    161            
     148            print 'General_mesh: Computing areas, normals and edgelengths'
     149
    162150        for i in range(N):
    163             if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N)
    164            
     151            if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N)
     152
    165153            x0, y0 = V[3*i, :]
    166154            x1, y1 = V[3*i+1, :]
    167             x2, y2 = V[3*i+2, :]           
     155            x2, y2 = V[3*i+2, :]
    168156
    169157            # 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]
     158            self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2
     159
     160            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2)
     161            msg += ' is degenerate:  area == %f' % self.areas[i]
    174162            assert self.areas[i] > 0.0, msg
    175 
    176163
    177164            # Normals
    … …  
    184171            #     the first vertex, etc)
    185172            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
    186 
    187             n0 = num.array([x2 - x1, y2 - y1], num.float)
     173            n0 = num.array([x2-x1, y2-y1], num.float)
    188174            l0 = num.sqrt(num.sum(n0**2))
    189175
    190             n1 = num.array([x0 - x2, y0 - y2], num.float)
     176            n1 = num.array([x0-x2, y0-y2], num.float)
    191177            l1 = num.sqrt(num.sum(n1**2))
    192178
    193             n2 = num.array([x1 - x0, y1 - y0], num.float)
     179            n2 = num.array([x1-x0, y1-y0], num.float)
    194180            l2 = num.sqrt(num.sum(n2**2))
    195181
    … …  
    207193            self.edgelengths[i, :] = [l0, l1, l2]
    208194
    209        
    210195        # Build structure listing which trianglse belong to which node.
    211         if verbose: print 'Building inverted triangle structure'         
     196        if verbose: print 'Building inverted triangle structure'
    212197        self.build_inverted_triangle_structure()
    213 
    214            
    215198
    216199    def __len__(self):
    217200        return self.number_of_triangles
    218    
    219201
    220202    def __repr__(self):
    221         return 'Mesh: %d vertices, %d triangles'\
    222                %(self.nodes.shape[0], len(self))
     203        return ('Mesh: %d vertices, %d triangles'
     204                % (self.nodes.shape[0], len(self)))
    223205
    224206    def get_normals(self):
    225207        """Return all normal vectors.
     208
    226209        Return normal vectors for all triangles as an Nx6 array
    227210        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
    228211        """
     212
    229213        return self.normals
    230 
    231214
    232215    def get_normal(self, i, j):
    233216        """Return normal vector j of the i'th triangle.
     217
    234218        Return value is the numeric array slice [x, y]
    235219        """
     220
    236221        return self.normals[i, 2*j:2*j+2]
    237222
    238223    def get_number_of_nodes(self):
    239224        return self.number_of_nodes
    240        
     225
    241226    def get_nodes(self, absolute=False):
    242227        """Return all nodes in mesh.
    … …  
    249234        are to be made absolute by taking georeference into account
    250235        Default is False as many parts of ANUGA expects relative coordinates.
    251         (To see which, switch to default absolute=True and run tests).       
     236        (To see which, switch to default absolute=True and run tests).
    252237        """
    253238
    … …  
    257242            if not self.geo_reference.is_absolute():
    258243                V = self.geo_reference.get_absolute(V)
    259                
     244
    260245        return V
    261    
    262     def get_node(self, i,
    263                  absolute=False):
     246
     247    def get_node(self, i, absolute=False):
    264248        """Return node coordinates for triangle i.
    265249
    … …  
    267251        are to be made absolute by taking georeference into account
    268252        Default is False as many parts of ANUGA expects relative coordinates.
    269         (To see which, switch to default absolute=True and run tests).       
    270         """
    271 
    272        
     253        (To see which, switch to default absolute=True and run tests).
     254        """
     255
    273256        V = self.nodes[i,:]
    274257        if absolute is True:
    275258            if not self.geo_reference.is_absolute():
    276                 return V + num.array([self.geo_reference.get_xllcorner(),
    277                                       self.geo_reference.get_yllcorner()], num.float)
    278             else:
    279                 return V
    280         else:
    281             return V
    282    
    283        
    284 
    285     def get_vertex_coordinates(self,
    286                                triangle_id=None,
    287                                absolute=False):
    288         """Return vertex coordinates for all triangles.
    289        
     259                V += num.array([self.geo_reference.get_xllcorner(),
     260                                self.geo_reference.get_yllcorner()], num.float)
     261        return V
     262
     263    def get_vertex_coordinates(self, triangle_id=None, absolute=False):
     264        """Return vertex coordinates for all triangles.
     265
    290266        Return all vertex coordinates for all triangles as a 3*M x 2 array
    291267        where the jth vertex of the ith triangle is located in row 3*i+j and
    … …  
    294270        if triangle_id is specified (an integer) the 3 vertex coordinates
    295271        for triangle_id are returned.
    296        
     272
    297273        Boolean keyword argument absolute determines whether coordinates
    298274        are to be made absolute by taking georeference into account
    299275        Default is False as many parts of ANUGA expects relative coordinates.
    300276        """
    301        
     277
    302278        V = self.vertex_coordinates
    303279
    304         if triangle_id is None:   
     280        if triangle_id is None:
    305281            if absolute is True:
    306282                if not self.geo_reference.is_absolute():
    307283                    V = self.geo_reference.get_absolute(V)
    308      
    309284            return V
    310285        else:
    … …  
    313288            assert int(i) == i, msg
    314289            assert 0 <= i < self.number_of_triangles
    315            
    316             i3 = 3*i 
     290
     291            i3 = 3*i
    317292            if absolute is True and not self.geo_reference.is_absolute():
    318293                offset=num.array([self.geo_reference.get_xllcorner(),
    … …  
    323298            else:
    324299                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
    325                
    326 
    327300
    328301    def get_vertex_coordinate(self, i, j, absolute=False):
    … …  
    333306        msg = 'vertex id j must be an integer in [0,1,2]'
    334307        assert j in [0,1,2], msg
    335        
    336         V = self.get_vertex_coordinates(triangle_id=i,
    337                                         absolute=absolute)
     308
     309        V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute)
    338310        return V[j,:]
    339    
    340 
    341311
    342312    def compute_vertex_coordinates(self):
    … …  
    358328        return vertex_coordinates
    359329
    360 
    361 
    362330    def get_triangles(self, indices=None):
    363331        """Get mesh triangles.
    … …  
    375343        if indices is None:
    376344            return self.triangles
    377             #indices = range(M)
    378345
    379346        return num.take(self.triangles, indices, axis=0)
    380    
    381 
    382347
    383348    def get_disconnected_triangles(self):
    … …  
    398363
    399364        The triangles created will have the format
    400 
    401365        [[0,1,2],
    402366         [3,4,5],
    403367         [6,7,8],
    404368         ...
    405          [3*M-3 3*M-2 3*M-1]]         
     369         [3*M-3 3*M-2 3*M-1]]
    406370        """
    407371
    408372        M = len(self) # Number of triangles
    409373        K = 3*M       # Total number of unique vertices
    410        
    411         T = num.reshape(num.arange(K, dtype=num.int), (M,3))
    412        
    413         return T     
    414 
    415    
    416 
    417     def get_unique_vertices(self,  indices=None):
     374        return num.reshape(num.arange(K, dtype=num.int), (M,3))
     375
     376    def get_unique_vertices(self, indices=None):
    418377        """FIXME(Ole): This function needs a docstring"""
    419378
    … …  
    421380        unique_verts = {}
    422381        for triangle in triangles:
    423             #print 'triangle(%s)=%s' % (type(triangle), str(triangle))
    424382            unique_verts[triangle[0]] = 0
    425383            unique_verts[triangle[1]] = 0
    … …  
    427385        return unique_verts.keys()
    428386
    429 
    430387    def get_triangles_and_vertices_per_node(self, node=None):
    431388        """Get triangles associated with given node.
    … …  
    441398            # Get index for this node
    442399            first = num.sum(self.number_of_triangles_per_node[:node])
    443            
     400
    444401            # Get number of triangles for this node
    445402            count = self.number_of_triangles_per_node[node]
    … …  
    459416            # working directly with the inverted triangle structure
    460417            for i in range(self.number_of_full_nodes):
    461                
    462418                L = self.get_triangles_and_vertices_per_node(node=i)
    463 
    464419                triangle_list.append(L)
    465420
    466421        return triangle_list
    467        
    468422
    469423    def build_inverted_triangle_structure(self):
    … …  
    475429        listing for each node how many triangles use it. N is the number of
    476430        nodes in mesh.
    477        
     431
    478432        vertex_value_indices: An array of length M listing indices into
    479433        triangles ordered by node number. The (triangle_id, vertex_id)
    … …  
    483437        used to average vertex values efficiently.
    484438
    485        
    486439        Example:
    487        
    488440        a = [0.0, 0.0] # node 0
    489441        b = [0.0, 2.0] # node 1
    … …  
    492444        e = [2.0, 2.0] # node 4
    493445        f = [4.0, 0.0] # node 5
    494 
    495446        nodes = array([a, b, c, d, e, f])
    496        
    497         #bac, bce, ecf, dbe, daf, dae
    498         triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])       
    499 
    500 
    501         For this structure
    502 
     447
     448        #                    bac,     bce,     ecf,     dbe
     449        triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
     450
     451        For this structure:
    503452        number_of_triangles_per_node = [1 3 3 1 3 1]
    504453        which means that node a has 1 triangle associated with it, node b
    505454        has 3, node has 3 and so on.
    506        
     455
    507456        vertex_value_indices = [ 1  0  3 10  2  4  7  9  5  6 11  8]
    508457        which reflects the fact that
    … …  
    518467                   and by triangle 2, vertex 0 (index = 6)
    519468                   and by triangle 3, vertex 2 (index = 11)
    520         node 5 is used by triangle 2, vertex 2 (index = 8)                   
    521        
     469        node 5 is used by triangle 2, vertex 2 (index = 8)
    522470
    523471        Preconditions:
    … …  
    526474        Postcondition:
    527475          self.number_of_triangles_per_node is built
    528           self.vertex_value_indices is built         
     476          self.vertex_value_indices is built
    529477        """
    530478
    … …  
    541489
    542490        # Register (triangle, vertex) indices for each node
    543         vertexlist = [None]*self.number_of_full_nodes
     491        vertexlist = [None] * self.number_of_full_nodes
    544492        for volume_id in range(self.number_of_full_triangles):
    545 
    546493            a = self.triangles[volume_id, 0]
    547494            b = self.triangles[volume_id, 1]
    548495            c = self.triangles[volume_id, 2]
    549496
    550             for vertex_id, node_id in enumerate([a,b,c]):
     497            for vertex_id, node_id in enumerate([a, b, c]):
    551498                if vertexlist[node_id] is None:
    552499                    vertexlist[node_id] = []
    553        
    554                 vertexlist[node_id].append( (volume_id, vertex_id) )
    555 
     500                vertexlist[node_id].append((volume_id, vertex_id))
    556501
    557502        # Build inverted triangle index array
    … …  
    561506                for volume_id, vertex_id in vertices:
    562507                    vertex_value_indices[k] = 3*volume_id + vertex_id
    563                                        
    564508                    k += 1
    565509
    … …  
    568512        self.vertex_value_indices = vertex_value_indices
    569513
    570 
    571        
    572 
    573514    def get_extent(self, absolute=False):
    574515        """Return min and max of all x and y coordinates
    … …  
    577518        are to be made absolute by taking georeference into account
    578519        """
    579        
    580 
    581520
    582521        C = self.get_vertex_coordinates(absolute=absolute)
    … …  
    588527        ymin = min(Y.flat)
    589528        ymax = max(Y.flat)
    590         #print "C",C
     529
    591530        return xmin, xmax, ymin, ymax
    592531
    593532    def get_areas(self):
    594         """Get areas of all individual triangles.
    595         """
    596         return self.areas       
     533        """Get areas of all individual triangles."""
     534
     535        return self.areas
    597536
    598537    def get_area(self):
    599         """Return total area of mesh
    600         """
     538        """Return total area of mesh"""
    601539
    602540        return num.sum(self.areas)
    … …  
    604542    def set_georeference(self, g):
    605543        self.geo_reference = g
    606        
     544
    607545    def get_georeference(self):
    608546        return self.geo_reference
    609        
    610        
    611                
     547

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

    r6304 r6428  
    1313
    1414import numpy as num
    15        
     15
    1616
    1717class Mesh(General_mesh):
    … …  
    9191                              verbose=verbose)
    9292
    93         if verbose: print 'Initialising mesh'         
     93        if verbose: print 'Initialising mesh'
    9494
    9595        N = len(self) #Number_of_triangles
    … …  
    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]
    … …  
    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):
    … …  
    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:
    … …  
    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
    … …  
    11411128            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
    11421129            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
    1143             d0 = num.sqrt(num.sum(z0**2)) 
     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:
    … …  
    11601147
    11611148
    1162             segment = ((x0,y0), (x1, y1))   
     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

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

    r6410 r6428  
    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
     59        #                        bac,     bce,     ecf,     dbe
    6160        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
    6261
    63         domain = General_mesh(nodes, triangles,
    64                        geo_reference = geo)
    65 
    66         verts = domain.get_vertex_coordinates(triangle_id=0)       
    67         msg = ("num.array([b,a,c])=\n%s\nshould be the same as 'verts'=\n%s"
     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"
    6866               % (str(num.array([b,a,c])), str(verts)))
    69         #self.assert_(num.allclose(num.array([b,a,c]), verts))
    70         assert num.allclose(num.array([b,a,c]), verts), msg
     67        self.failUnless(num.allclose(num.array([b,a,c]), verts), msg)
    7168
    7269        verts = domain.get_vertex_coordinates(triangle_id=0)       
    73         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
    7485        verts = domain.get_vertex_coordinates(triangle_id=0,
    7586                                              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)))
    7692        self.assert_(num.allclose(num.array([nodes_absolute[1],
    77                                      nodes_absolute[0],
    78                                      nodes_absolute[2]]), verts))
    79         verts = domain.get_vertex_coordinates(triangle_id=0,
    80                                               absolute=True)       
    81         self.assert_(num.allclose(num.array([nodes_absolute[1],
    82                                      nodes_absolute[0],
    83                                      nodes_absolute[2]]), verts))
    84        
    85        
     93                                             nodes_absolute[0],
     94                                             nodes_absolute[2]]),
     95                                  verts), msg)
    8696
    8797    def test_get_vertex_coordinates_triangle_id(self):
    … …  
    283293        nodes_absolute = geo.get_absolute(nodes)
    284294       
    285         #bac, bce, ecf, dbe, daf, dae
     295        #                        bac,     bce,     ecf,     dbe
    286296        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
    287297
    288         domain = General_mesh(nodes, triangles,
    289                        geo_reference = geo)
     298        domain = General_mesh(nodes, triangles, geo_reference = geo)
    290299        node = domain.get_node(2)       
    291         #self.assertEqual(c, node)
    292         msg = ('\nc=%s\nmode=%s' % (str(c), str(node)))
     300        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
    293301        self.failUnless(num.alltrue(c == node), msg)
     302
     303        # repeat get_node(), see if result same
     304        node = domain.get_node(2)       
     305        msg = ('\nc=%s\nnode=%s' % (str(c), str(node)))
     306        self.failUnless(num.alltrue(c == node), msg)
    294307       
    295308        node = domain.get_node(2, absolute=True)     
    296         #self.assertEqual(nodes_absolute[2], node)
    297         self.failUnless(num.alltrue(nodes_absolute[2] == node))
    298        
     309        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
     310               % (str(nodes_absolute[2]), str(node)))
     311        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
     312       
     313        # repeat get_node(absolute=True), see if result same
    299314        node = domain.get_node(2, absolute=True)     
    300         #self.assertEqual(nodes_absolute[2], node)
    301         self.failUnless(num.alltrue(nodes_absolute[2] == node))
     315        msg = ('\nnodes_absolute[2]=%s\nnode=%s'
     316               % (str(nodes_absolute[2]), str(node)))
     317        self.failUnless(num.alltrue(nodes_absolute[2] == node), msg)
    302318       
    303319
    … …  
    337353        self.failUnlessRaises(AssertionError, General_mesh,
    338354                              nodes, triangles, geo_reference=geo)
     355
     356    def test_raw_change_points_geo_ref(self):
     357        x0 = 1000.0
     358        y0 = 2000.0
     359        geo = Geo_reference(56, x0, y0)
     360       
    339361       
    340362
    … …  
    343365
    344366if __name__ == "__main__":
    345     suite = unittest.makeSuite(Test_General_Mesh,'test')
     367    suite = unittest.makeSuite(Test_General_Mesh, 'test')
     368    #suite = unittest.makeSuite(Test_General_Mesh, 'test_get_vertex_coordinates_with_geo_ref')
    346369    runner = unittest.TextTestRunner()
    347370    runner.run(suite)

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

    r6410 r6428  
    11691169        """
    11701170       
    1171        
    11721171        # test
    11731172        a = [0.0, 0.0]
    … …  
    11761175
    11771176        absolute_points = [a, b, c]
    1178         vertices = [[0,1,2]]
    1179        
    1180         geo = Geo_reference(56,67,-56)
     1177        vertices = [[0, 1, 2]]
     1178       
     1179        geo = Geo_reference(56, 67, -56)
    11811180
    11821181        relative_points = geo.change_points_geo_ref(absolute_points)
    1183 
    1184         #print 'Relative', relative_points
    1185         #print 'Absolute', absolute_points       
    11861182
    11871183        mesh = Mesh(relative_points, vertices, geo_reference=geo)
    … …  
    18421838
    18431839if __name__ == "__main__":
    1844     suite = unittest.makeSuite(Test_Mesh,'test')
     1840    suite = unittest.makeSuite(Test_Mesh, 'test')
    18451841    runner = unittest.TextTestRunner()
    18461842    runner.run(suite)

  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r6360 r6428  
    101101            self.read_ASCII(ASCIIFile, read_title=read_title)
    102102
    103 #    # Might be better to have this method instead of the following 3.
    104 #    def get_origin(self):
    105 #        return (self.zone, self.xllcorner, self.yllcorner)
    106 
    107103    ##
    108104    # @brief Get the X coordinate of the origin of this georef.
    … …  
    238234################################################################################
    239235
     236    ##
     237    # @brief Change points to be absolute wrt new georef 'points_geo_ref'.
     238    # @param points The points to change.
     239    # @param points_geo_ref The new georef to make points absolute wrt.
     240    # @return The changed points.
     241    # @note If 'points' is a list then a changed list is returned.
     242    # @note The input points data is changed.
    240243    def change_points_geo_ref(self, points, points_geo_ref=None):
     244        """Change the geo reference of a list or numeric array of points to
     245        be this reference.(The reference used for this object)
     246        If the points do not have a geo ref, assume 'absolute' input values
    241247        """
    242         Change the geo reference of a list or numeric array of points to
    243         be this reference.(The reference used for this object)
    244         If the points do not have a geo ref, assume 'absolute' values
    245         """
    246 
    247         is_list = False
    248         if type(points) == types.ListType:
    249             is_list = True
     248
     249        # remember if we got a list
     250        is_list = isinstance(points, list)
    250251
    251252        points = ensure_numeric(points, num.float)
    252253
     254        # sanity checks
    253255        if len(points.shape) == 1:
    254256            # One point has been passed
    … …  
    288290    # @return True if ???
    289291    def is_absolute(self):
    290         """Return True if xllcorner==yllcorner==0
    291         indicating that points in question are absolute.
    292         """
     292        """Return True if xllcorner==yllcorner==0"""
    293293
    294294        return num.allclose([self.xllcorner, self.yllcorner], 0)
    … …  
    298298    # @param points
    299299    # @return
    300     # @note
     300    # @note This method changes the input points!
    301301    def get_absolute(self, points):
    302         """Given a set of points geo referenced to this instance,
     302        """Given a set of points geo referenced to this instance
     303
    303304        return the points as absolute values.
    304305        """
    305306
    306         #if self.is_absolute:
    307         #    return points
    308 
    309307        # remember if we got a list
    310         is_list = False
    311         if type(points) == types.ListType:
    312             is_list = True
    313 
     308        is_list = isinstance(points, list)
     309
     310        # convert to numeric array
    314311        points = ensure_numeric(points, num.float)
     312
     313        # sanity checks
    315314        if len(points.shape) == 1:
    316315            #One point has been passed
    … …  
    326325
    327326        # Add geo ref to points
    328         #if not self.is_absolute:
    329327        if not self.is_absolute():
    330328            points[:,0] += self.xllcorner

  • branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

    r6410 r6428  
    503503                        'bad shape did not raise error!')
    504504            os.remove(point_file)
     505
     506    def test_functionality_get_absolute(self):
     507        x0 = 1000.0
     508        y0 = 2000.0
     509        geo = Geo_reference(56, x0, y0)
     510
     511        # iterable points (*not* num.array())
     512        points = ((2,3), (3,1), (5,2))
     513        abs_points = geo.get_absolute(points)
     514        # check we haven't changed 'points' itself
     515        self.failIf(num.alltrue(abs_points == points))
     516        new_points = abs_points.copy()
     517        new_points[:,0] -= x0
     518        new_points[:,1] -= y0
     519        self.failUnless(num.alltrue(new_points == points))
     520
     521        # points in num.array()
     522        points = num.array(((2,3), (3,1), (5,2)), num.float)
     523        abs_points = geo.get_absolute(points)
     524        # check we haven't changed 'points' itself
     525        self.failIf(num.alltrue(abs_points == points))
     526        new_points = abs_points.copy()
     527        new_points[:,0] -= x0
     528        new_points[:,1] -= y0
     529        self.failUnless(num.alltrue(new_points == points))
     530
    505531       
    506532#-------------------------------------------------------------
    507533
    508534if __name__ == "__main__":
    509     suite = unittest.makeSuite(geo_referenceTestCase,'test')
     535    suite = unittest.makeSuite(geo_referenceTestCase, 'test')
     536    #suite = unittest.makeSuite(geo_referenceTestCase, 'test_functionality_get_absolute')
    510537    runner = unittest.TextTestRunner() #verbosity=2)
    511538    runner.run(suite)

  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6410 r6428  
    305305                    t = self.interpolate_block(f, point_coordinates[start:end],
    306306                                               verbose=verbose)
    307                     z = num.concatenate((z, t))
     307                    z = num.concatenate((z, t), axis=0)    #??default#
    308308                    start = end
    309309
    … …  
    311311                t = self.interpolate_block(f, point_coordinates[start:end],
    312312                                           verbose=verbose)
    313                 z = num.concatenate((z, t))
     313                z = num.concatenate((z, t), axis=0)    #??default#
    314314        return z
    315315
    … …  
    11891189
    11901190    #Add the x and y together
    1191     vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),axis=1)
     1191    vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),
     1192                                         axis=1)
    11921193
    11931194    #Will return the quantity values at the specified times and locations

  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6410 r6428  
    1 """Class Geospatial_data - Manipulation of locations on the planet and
    2 associated attributes.
    3 
     1"""Class Geospatial_data
     2
     3Manipulation of locations on the planet and associated attributes.
    44"""
    55
    … …  
    1111from RandomArray import randint, seed, get_seed
    1212from copy import deepcopy
     13
    1314from Scientific.IO.NetCDF import NetCDFFile
    14 
    1515import numpy as num
    1616
    … …  
    2525from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    2626from anuga.config import netcdf_float
     27
    2728
    2829DEFAULT_ATTRIBUTE = 'elevation'
    … …  
    139140
    140141        self.set_verbose(verbose)
    141         self.geo_reference = None # create the attribute
     142        self.geo_reference = None
    142143        self.file_name = file_name
    143144
    … …  
    148149
    149150        if file_name is None:
    150             if latitudes is not None \
    151                or longitudes is not None \
    152                or points_are_lats_longs:
     151            if (latitudes is not None or longitudes is not None
     152                    or points_are_lats_longs):
    153153                data_points, geo_reference = \
    154154                    _set_using_lat_long(latitudes=latitudes,
    … …  
    156156                                        geo_reference=geo_reference,
    157157                                        data_points=data_points,
    158                                         points_are_lats_longs=\
    159                                         points_are_lats_longs)
     158                                        points_are_lats_longs=
     159                                            points_are_lats_longs)
    160160            self.check_data_points(data_points)
    161161            self.set_attributes(attributes)
    … …  
    213213            msg = 'There is no data or file provided!'
    214214            raise ValueError, msg
    215 
    216215        else:
    217216            self.data_points = ensure_numeric(data_points)
    … …  
    257256        """
    258257
    259         from anuga.coordinate_transforms.geo_reference import Geo_reference
    260 
    261258        if geo_reference is None:
    262259            # Use default - points are in absolute coordinates
    … …  
    269266            # FIXME (Ole): This exception will be raised even
    270267            # if geo_reference is None. Is that the intent Duncan?
    271             msg = 'Argument geo_reference must be a valid Geo_reference '
    272             msg += 'object or None.'
     268            msg = ('Argument geo_reference must be a valid Geo_reference '
     269                   'object or None.')
    273270            raise Expection, msg
    274271
    … …  
    378375    # @param isSouthHemisphere If True, return lat/lon points in S.Hemi.
    379376    # @return A set of data points, in appropriate form.
    380     def get_data_points(self, absolute=True, geo_reference=None,
    381                         as_lat_long=False, isSouthHemisphere=True):
     377    def get_data_points(self,
     378                        absolute=True,
     379                        geo_reference=None,
     380                        as_lat_long=False,
     381                        isSouthHemisphere=True):
    382382        """Get coordinates for all data points as an Nx2 array
    383383
    … …  
    460460    # @return The new geospatial object.
    461461    def __add__(self, other):
    462         """Returns the addition of 2 geospatical objects,
     462        """Returns the addition of 2 geospatial objects,
    463463        objects are concatencated to the end of each other
    464464
    … …  
    488488            if self.attributes is None:
    489489                if other.attributes is not None:
    490                     msg = 'Geospatial data must have the same '
    491                     msg += 'attributes to allow addition.'
     490                    msg = ('Geospatial data must have the same '
     491                           'attributes to allow addition.')
    492492                    raise Exception, msg
    493493
    … …  
    499499                        attrib1 = self.attributes[x]
    500500                        attrib2 = other.attributes[x]
    501                         new_attributes[x] = num.concatenate((attrib1, attrib2))
     501                        new_attributes[x] = num.concatenate((attrib1, attrib2),
     502                                                            axis=0) #??default#
    502503                    else:
    503                         msg = 'Geospatial data must have the same \n'
    504                         msg += 'attributes to allow addition.'
     504                        msg = ('Geospatial data must have the same '
     505                               'attributes to allow addition.')
    505506                        raise Exception, msg
    506507        else:
    … …  
    523524        return self + other
    524525
    525     ###
    526     #  IMPORT/EXPORT POINTS FILES
    527     ###
     526################################################################################
     527#  IMPORT/EXPORT POINTS FILES
     528################################################################################
    528529
    529530    ##
    … …  
    560561            except IOError, e:
    561562                # This should only be if a file is not found
    562                 msg = 'Could not open file %s. ' % file_name
    563                 msg += 'Check the file location.'
     563                msg = ('Could not open file %s. Check the file location.'
     564                       % file_name)
    564565                raise IOError, msg
    565566            except SyntaxError, e:
    566567                # This should only be if there is a format error
    567                 msg = 'Problem with format of file %s. \n' %file_name
    568                 msg += Error_message['IOError']
     568                msg = ('Problem with format of file %s.\n%s'
     569                       % (file_name, Error_message['IOError']))
    569570                raise SyntaxError, msg
    570571        else:
    571             msg = 'Extension %s is unknown' %file_name[-4:]
     572            msg = 'Extension %s is unknown' % file_name[-4:]
    572573            raise IOError, msg
    573574
    … …  
    614615                                self.get_all_attributes(),
    615616                                self.get_geo_reference())
    616 
    617617        elif file_name[-4:] == ".txt" or file_name[-4:] == ".csv":
    618618            msg = "ERROR: trying to write a .txt file with relative data."
    … …  
    621621                            self.get_data_points(absolute=True,
    622622                                                 as_lat_long=as_lat_long,
    623                                   isSouthHemisphere=isSouthHemisphere),
     623                                           isSouthHemisphere=isSouthHemisphere),
    624624                            self.get_all_attributes(),
    625625                            as_lat_long=as_lat_long)
    626 
    627626        elif file_name[-4:] == ".urs" :
    628627            msg = "ERROR: Can not write a .urs file as a relative file."
    … …  
    630629            _write_urs_file(file_name,
    631630                            self.get_data_points(as_lat_long=True,
    632                                   isSouthHemisphere=isSouthHemisphere))
    633 
     631                                           isSouthHemisphere=isSouthHemisphere))
    634632        else:
    635633            msg = 'Unknown file type %s ' %file_name
    … …  
    694692        random_list = []
    695693        remainder_list = []
    696         new_size = round(factor*self_size)
     694        new_size = round(factor * self_size)
    697695
    698696        # Find unique random numbers
    … …  
    713711        # Set seed if provided, mainly important for unit test!
    714712        # plus recalcule seed when no seed provided.
    715         if seed_num != None:
     713        if seed_num is not None:
    716714            seed(seed_num, seed_num)
    717715        else:
    … …  
    734732
    735733        # pops array index (random_num) from remainder_list
    736         # (which starts as the
    737         # total_list and appends to random_list
     734        # (which starts as the total_list and appends to random_list)
    738735        random_num_len = len(random_num)
    739736        for i in random_num:
    … …  
    746743
    747744        # FIXME: move to tests, it might take a long time
    748         # then create an array of random lenght between 500 and 1000,
     745        # then create an array of random length between 500 and 1000,
    749746        # and use a random factor between 0 and 1
    750747        # setup for assertion
    … …  
    752749        test_total.extend(remainder_list)
    753750        test_total.sort()
    754         msg = 'The two random lists made from the original list when added ' \
    755               'together DO NOT equal the original list'
     751        msg = ('The two random lists made from the original list when added '
     752               'together DO NOT equal the original list')
    756753        assert total_list == test_total, msg
    757754
    … …  
    770767        file pointer position
    771768        """
    772         from Scientific.IO.NetCDF import NetCDFFile
     769
    773770        # FIXME - what to do if the file isn't there
    774771
    … …  
    799796            self.number_of_blocks = self.number_of_points/self.max_read_lines
    800797            # This computes the number of full blocks. The last block may be
    801             # smaller and won't be ircluded in this estimate.
     798            # smaller and won't be included in this estimate.
    802799
    803800            if self.verbose is True:
    804                 print 'Reading %d points (in ~%d blocks) from file %s. ' \
    805                       % (self.number_of_points,
    806                          self.number_of_blocks,
    807                          self.file_name),
    808                 print 'Each block consists of %d data points' \
    809                       % self.max_read_lines
    810 
     801                print ('Reading %d points (in ~%d blocks) from file %s. '
     802                       % (self.number_of_points, self.number_of_blocks,
     803                          self.file_name)),
     804                print ('Each block consists of %d data points'
     805                       % self.max_read_lines)
    811806        else:
    812807            # Assume the file is a csv file
    … …  
    839834
    840835            if self.verbose is True:
    841                 if self.show_verbose >= self.start_row \
    842                    and self.show_verbose < fin_row:
    843                     print 'Reading block %d (points %d to %d) out of %d'\
    844                           %(self.block_number,
    845                             self.start_row,
    846                             fin_row,
    847                             self.number_of_blocks)
     836                if (self.show_verbose >= self.start_row
     837                    and self.show_verbose < fin_row):
     838                    print ('Reading block %d (points %d to %d) out of %d'
     839                           % (self.block_number, self.start_row,
     840                              fin_row, self.number_of_blocks))
    848841
    849842                    self.show_verbose += max(self.max_read_lines,
    850843                                             self.verbose_block_size)
    851 
    852844
    853845            # Read next block
    … …  
    861853
    862854            self.block_number += 1
    863 
    864855        else:
    865856            # Assume the file is a csv file
    … …  
    868859                 att_dict,
    869860                 geo_ref,
    870                  self.file_pointer) = \
    871                         _read_csv_file_blocking(self.file_pointer,
    872                                                 self.header[:],
    873                                                 max_read_lines=\
    874                                                     self.max_read_lines,
    875                                                 verbose=self.verbose)
     861                 self.file_pointer) = _read_csv_file_blocking(self.file_pointer,
     862                                                              self.header[:],
     863                                                              max_read_lines= \
     864                                                           self.max_read_lines,
     865                                                              verbose= \
     866                                                                  self.verbose)
    876867
    877868                # Check that the zones haven't changed.
    … …  
    880871                    self.blocking_georef = geo_ref
    881872                elif self.blocking_georef is not None:
    882                     msg = 'Geo reference given, then not given.'
    883                     msg += ' This should not happen.'
     873                    msg = ('Geo reference given, then not given.'
     874                           ' This should not happen.')
    884875                    raise ValueError, msg
    885876                geo = Geospatial_data(pointlist, att_dict, geo_ref)
    … …  
    899890                del self.file_pointer
    900891                # This should only be if there is a format error
    901                 msg = 'Could not open file %s. \n' % self.file_name
    902                 msg += Error_message['IOError']
     892                msg = ('Could not open file %s.\n%s'
     893                       % (self.file_name, Error_message['IOError']))
    903894                raise SyntaxError, msg
    904895        return geo
    905 
    906896
    907897##################### Error messages ###########
    908898Error_message = {}
    909899Em = Error_message
    910 Em['IOError'] = "NOTE: The format for a comma separated .txt/.csv file is:\n"
    911 Em['IOError'] += "        1st line:     [column names]\n"
    912 Em['IOError'] += "        other lines:  [x value], [y value], [attributes]\n"
    913 Em['IOError'] += "\n"
    914 Em['IOError'] += "           for example:\n"
    915 Em['IOError'] += "           x, y, elevation, friction\n"
    916 Em['IOError'] += "           0.6, 0.7, 4.9, 0.3\n"
    917 Em['IOError'] += "           1.9, 2.8, 5, 0.3\n"
    918 Em['IOError'] += "           2.7, 2.4, 5.2, 0.3\n"
    919 Em['IOError'] += "\n"
    920 Em['IOError'] += "The first two columns are assumed to be x, y coordinates.\n"
    921 Em['IOError'] += "The attribute values must be numeric.\n"
     900Em['IOError'] = ('NOTE: The format for a comma separated .txt/.csv file is:\n'
     901                 '        1st line:     [column names]\n'
     902                 '        other lines:  [x value], [y value], [attributes]\n'
     903                 '\n'
     904                 '           for example:\n'
     905                 '           x, y, elevation, friction\n'
     906                 '           0.6, 0.7, 4.9, 0.3\n'
     907                 '           1.9, 2.8, 5, 0.3\n'
     908                 '           2.7, 2.4, 5.2, 0.3\n'
     909                 '\n'
     910                 'The first two columns are assumed to be x, y coordinates.\n'
     911                 'The attribute values must be numeric.\n')
    922912
    923913##
    … …  
    936926
    937927    if geo_reference is not None:
    938         msg = "A georeference is specified yet latitude and longitude " \
    939               "are also specified!"
     928        msg = ('A georeference is specified yet latitude and longitude '
     929               'are also specified!')
    940930        raise ValueError, msg
    941931
    942932    if data_points is not None and not points_are_lats_longs:
    943         msg = "Data points are specified yet latitude and longitude are " \
    944               "also specified."
     933        msg = ('Data points are specified yet latitude and longitude are '
     934               'also specified.')
    945935        raise ValueError, msg
    946936
    … …  
    983973    dic['attributelist']['elevation'] = [[7.0,5.0]]
    984974    """
    985 
    986     from Scientific.IO.NetCDF import NetCDFFile
    987975
    988976    if verbose: print 'Reading ', file_name
    … …  
    11271115            if len(header) != len(numbers):
    11281116                file_pointer.close()
    1129                 msg = "File load error. " \
    1130                       "There might be a problem with the file header."
     1117                msg = ('File load error. '
     1118                       'There might be a problem with the file header.')
    11311119                raise SyntaxError, msg
    11321120            for i,n in enumerate(numbers):
    11331121                n.strip()
    11341122                if n != '\n' and n != '':
    1135                     #attributes.append(float(n))
    11361123                    att_dict.setdefault(header[i],[]).append(float(n))
    11371124        except ValueError:
    … …  
    11741161# @note Will throw IOError and AttributeError exceptions.
    11751162def _read_pts_file_header(fid, verbose=False):
    1176     """Read the geo_reference and number_of_points from a .pts file"""
     1163    '''Read the geo_reference and number_of_points from a .pts file'''
    11771164
    11781165    keys = fid.variables.keys()
    … …  
    12021189# @return Tuple of (pointlist, attributes).
    12031190def _read_pts_file_blocking(fid, start_row, fin_row, keys):
    1204     """Read the body of a .csv file."""
     1191    '''Read the body of a .csv file.'''
    12051192
    12061193    pointlist = num.array(fid.variables['points'][start_row:fin_row])
    … …  
    12401227    """
    12411228
    1242     from Scientific.IO.NetCDF import NetCDFFile
    1243 
    12441229    # NetCDF file definition
    12451230    outfile = NetCDFFile(file_name, netcdf_mode_w)
    … …  
    12561241
    12571242    # Variable definition
    1258     outfile.createVariable('points', netcdf_float, ('number_of_points',
    1259                                                    'number_of_dimensions'))
     1243    outfile.createVariable('points', netcdf_float,
     1244                           ('number_of_points', 'number_of_dimensions'))
    12601245
    12611246    # create variables
    … …  
    13851370    """Convert points_dictionary to geospatial data object"""
    13861371
    1387     msg = 'Points dictionary must have key pointlist'
     1372    msg = "Points dictionary must have key 'pointlist'"
    13881373    assert points_dictionary.has_key('pointlist'), msg
    13891374
    1390     msg = 'Points dictionary must have key attributelist'
     1375    msg = "Points dictionary must have key 'attributelist'"
    13911376    assert points_dictionary.has_key('attributelist'), msg
    13921377
    … …  
    13981383    return Geospatial_data(points_dictionary['pointlist'],
    13991384                           points_dictionary['attributelist'],
    1400                            geo_reference = geo)
     1385                           geo_reference=geo)
    14011386
    14021387
    … …  
    14361421    # Sort of geo_reference and convert points
    14371422    if geo_reference is None:
    1438         geo = None # Geo_reference()
     1423        geo = None    # Geo_reference()
    14391424    else:
    14401425        if isinstance(geo_reference, Geo_reference):
    … …  
    15721557    from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    15731558
    1574 
    15751559    attribute_smoothed = 'elevation'
    15761560
    … …  
    15791563        mesh_file = 'temp.msh'
    15801564
    1581         if north_boundary is None or south_boundary is None \
    1582            or east_boundary is None or west_boundary is None:
     1565        if (north_boundary is None or south_boundary is None
     1566            or east_boundary is None or west_boundary is None):
    15831567            no_boundary = True
    15841568        else:
    … …  
    15891573            raise Expection, msg
    15901574
    1591         poly_topo = [[east_boundary,south_boundary],
    1592                      [east_boundary,north_boundary],
    1593                      [west_boundary,north_boundary],
    1594                      [west_boundary,south_boundary]]
     1575        poly_topo = [[east_boundary, south_boundary],
     1576                     [east_boundary, north_boundary],
     1577                     [west_boundary, north_boundary],
     1578                     [west_boundary, south_boundary]]
    15951579
    15961580        create_mesh_from_regions(poly_topo,
    … …  
    16221606    if alpha_list == None:
    16231607        alphas = [0.001,0.01,100]
    1624         #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
    1625          #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     1608        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
     1609        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
    16261610    else:
    16271611        alphas = alpha_list
    … …  
    16471631        # add G_other data to domains with different alphas
    16481632        if verbose:
    1649             print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
     1633            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
    16501634        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
    16511635        if verbose: print domain.statistics()
    … …  
    16711655        sample_cov = cov(elevation_sample)
    16721656        ele_cov = cov(elevation_sample - elevation_predicted)
    1673         normal_cov[i,:] = [alpha,ele_cov / sample_cov]
     1657        normal_cov[i,:] = [alpha, ele_cov / sample_cov]
    16741658
    16751659        if verbose:
    … …  
    16941678        print 'Final results:'
    16951679        for i, alpha in enumerate(alphas):
    1696             print 'covariance for alpha %s = %s ' \
    1697                   % (normal_cov[i][0], normal_cov[i][1])
    1698         print '\n Optimal alpha is: %s ' \
    1699               % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0]
     1680            print ('covariance for alpha %s = %s '
     1681                   % (normal_cov[i][0], normal_cov[i][1]))
     1682        print ('\nOptimal alpha is: %s '
     1683               % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0])
    17001684
    17011685    # covariance and optimal alpha
    … …  
    17791763    from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    17801764
    1781 
    17821765    attribute_smoothed = 'elevation'
    17831766
    … …  
    17851768        mesh_file = 'temp.msh'
    17861769
    1787         if north_boundary is None or south_boundary is None \
    1788            or east_boundary is None or west_boundary is None:
     1770        if (north_boundary is None or south_boundary is None
     1771            or east_boundary is None or west_boundary is None):
    17891772            no_boundary = True
    17901773        else:
    … …  
    17951778            raise Expection, msg
    17961779
    1797         poly_topo = [[east_boundary,south_boundary],
    1798                      [east_boundary,north_boundary],
    1799                      [west_boundary,north_boundary],
    1800                      [west_boundary,south_boundary]]
     1780        poly_topo = [[east_boundary, south_boundary],
     1781                     [east_boundary, north_boundary],
     1782                     [west_boundary, north_boundary],
     1783                     [west_boundary, south_boundary]]
    18011784
    18021785        create_mesh_from_regions(poly_topo,
    … …  
    18221805    points = G_small.get_data_points()
    18231806
    1824     # FIXME: Remove points outside boundary polygon
    1825 #    print 'new point',len(points)
    1826 #
    1827 #    new_points=[]
    1828 #    new_points=array([],dtype=float)
    1829 #    new_points=resize(new_points,(len(points),2))
    1830 #    print "BOUNDARY", boundary_poly
    1831 #    for i,point in enumerate(points):
    1832 #        if is_inside_polygon(point,boundary_poly, verbose=True):
    1833 #            new_points[i] = point
    1834 #            print"WOW",i,new_points[i]
    1835 #    points = new_points
    1836 
    18371807    if verbose: print "Number of points in sample to compare: ", len(points)
    18381808
    18391809    if alpha_list == None:
    18401810        alphas = [0.001,0.01,100]
    1841         #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
    1842          #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     1811        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
     1812        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
    18431813    else:
    18441814        alphas = alpha_list
    … …  
    18511821        # add G_other data to domains with different alphas
    18521822        if verbose:
    1853             print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
     1823            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
    18541824        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
    18551825        if verbose: print domain.statistics()
    … …  
    18781848    if verbose:
    18791849        print 'Determine difference between predicted results and actual data'
     1850
    18801851    for i, alpha in enumerate(domains):
    18811852        if verbose: print'Alpha =', alpha

  • branches/numpy/anuga/geospatial_data/test_geospatial_data.py

    r6410 r6428  
    341341        P = G.get_data_points(absolute=True)
    342342
    343         assert num.allclose(P, num.concatenate( (P1,P2) ))
     343        assert num.allclose(P, num.concatenate((P1,P2), axis=0))    #??default#
    344344        msg = ('P=\n%s\nshould be close to\n'
    345345               '[[2.0, 4.1], [4.0, 7.3],\n'
    … …  
    389389        P = G.get_data_points(absolute=True)
    390390
    391         assert num.allclose(P, num.concatenate((points1, points2)))
     391        assert num.allclose(P, num.concatenate((points1, points2), axis=0)) #??default#
    392392
    393393    def test_add_with_None(self):

  • branches/numpy/anuga/load_mesh/loadASCII.py

    r6410 r6428  
    6666from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    6767from anuga.config import netcdf_float, netcdf_char, netcdf_int
     68from anuga.utilities.system_tools import *
    6869
    6970from Scientific.IO.NetCDF import NetCDFFile
    … …  
    624625            num.array(mesh['vertex_attributes'], num.float)
    625626    mesh['vertex_attribute_titles'] = \
    626         num.array(mesh['vertex_attribute_titles'], num.character)
     627        num.array(string_to_char(mesh['vertex_attribute_titles']), num.character)
    627628    mesh['segments'] = num.array(mesh['segments'], IntType)
    628     mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character)
     629    mesh['segment_tags'] = num.array(string_to_char(mesh['segment_tags']),
     630                                     num.character)
    629631    mesh['triangles'] = num.array(mesh['triangles'], IntType)
    630     mesh['triangle_tags'] = num.array(mesh['triangle_tags'])
     632    mesh['triangle_tags'] = num.array(string_to_char(mesh['triangle_tags']),
     633                                      num.character)
    631634    mesh['triangle_neighbors'] = \
    632635        num.array(mesh['triangle_neighbors'], IntType)
    … …  
    637640    mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType)
    638641    mesh['outline_segment_tags'] = \
    639         num.array(mesh['outline_segment_tags'], num.character)
     642        num.array(string_to_char(mesh['outline_segment_tags']), num.character)
    640643    mesh['holes'] = num.array(mesh['holes'], num.float)
    641644    mesh['regions'] = num.array(mesh['regions'], num.float)
    642     mesh['region_tags'] = num.array(mesh['region_tags'], num.character)
     645    mesh['region_tags'] = num.array(string_to_char(mesh['region_tags']), num.character)
    643646    mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float)
    644647
    … …  
    697700                               ('num_of_segments', 'num_of_segment_ends'))
    698701        outfile.variables['segments'][:] = mesh['segments']
    699         if mesh['segment_tags'].shape[1] > 0:
     702        if mesh['segment_tags'].shape[0] > 0:
    700703            outfile.createDimension('num_of_segment_tag_chars',
    701704                                    mesh['segment_tags'].shape[1])
    … …  
    747750                                'num_of_segment_ends'))
    748751        outfile.variables['outline_segments'][:] = mesh['outline_segments']
    749         if (mesh['outline_segment_tags'].shape[1] > 0):
     752        if mesh['outline_segment_tags'].shape[1] > 0:
    750753            outfile.createDimension('num_of_outline_segment_tag_chars',
    751754                                    mesh['outline_segment_tags'].shape[1])
    … …  
    817820    try:
    818821        titles = fid.variables['vertex_attribute_titles'][:]
    819         for i, title in enumerate(titles):
    820             mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
    821     except KeyError:
     822        mesh['vertex_attribute_titles'] = [x.tostring().strip() for x in titles]
     823    except:
    822824        pass
    823825
    … …  
    827829        mesh['segments'] = num.array([], num.int)      #array default#
    828830
    829     mesh['segment_tags'] =[]
     831    mesh['segment_tags'] = []
    830832    try:
    831833        tags = fid.variables['segment_tags'][:]
    832         for i, tag in enumerate(tags):
    833             mesh['segment_tags'].append(tags[i].tostring().strip())
    834     except KeyError:
    835         for ob in mesh['segments']:
    836             mesh['segment_tags'].append('')
     834        mesh['segment_tags'] = [x.tostring().strip() for x in tags]
     835    except:
     836        pass
    837837
    838838    try:
    … …  
    846846    try:
    847847        tags = fid.variables['triangle_tags'][:]
    848         for i, tag in enumerate(tags):
    849             mesh['triangle_tags'].append(tags[i].tostring().strip())
     848        mesh['triangle_tags'] = [x.tostring().strip() for x in tags]
    850849    except KeyError:
    851850        for ob in mesh['triangles']:

  • branches/numpy/anuga/load_mesh/test_loadASCII.py

    r6410 r6428  
    209209        os.remove(fileName)
    210210        dict = self.dict
    211         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file')
     211        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_tsh_file')
    212212
    213213    def test_read_write_tsh_fileII(self):
    … …  
    227227        os.remove(fileName)
    228228        dict = self.blank_dict
    229         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_fileIII')
     229        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_tsh_fileIII')
    230230
    231231    def test_read_write_tsh_file4(self):
    … …  
    236236        os.remove(fileName)
    237237        dict = self.seg_dict
    238         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file4')
     238        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_tsh_file4')
    239239
    240240    def test_read_write_tsh_file5(self):
    … …  
    244244        loaded_dict = import_mesh_file(fileName)
    245245        dict = self.triangle_tags_dict
    246         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file5')
     246        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_tsh_file5')
    247247        os.remove(fileName)
    248248
    … …  
    253253        loaded_dict = import_mesh_file(fileName)
    254254        dict = self.tri_dict
    255         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file6')
     255        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_tsh_file6')
    256256        os.remove(fileName)
    257257
    … …  
    312312  ############### .MSH ##########
    313313
    314     def test_read_write_msh_file(self):
     314    def test_read_write_msh_file1(self):
    315315        dict = self.dict.copy()
    316316        fileName = tempfile.mktemp('.msh')
    … …  
    319319        os.remove(fileName)
    320320        dict = self.dict
    321         self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file')
     321        self.check_mesh_dicts(loaded_dict, dict, 'test_read_write_msh_file1')
    322322
    323323    def test_read_write_msh_fileII(self):
    … …  
    399399                            num.array(loaded_dict['segments']))
    400400
     401        self.failUnlessEqual(dict['triangle_tags'], loaded_dict['triangle_tags'])
     402
    401403        for ob, ldob in map(None, dict['triangle_tags'],
    402404                            loaded_dict['triangle_tags']):
     405            msg = ('ob=\n%s\nshould be same as ldob=\n%s' % (str(ob), str(ldob)))
    403406            self.failUnless(ob == ldob,
    404                             fail_string + ' failed!! Test triangle_tags')
     407                            fail_string + ' failed\n' + msg)
    405408
    406409        # A bit hacky
    407         self.failUnless((loaded_dict['vertex_attributes'] ==
    408                          dict['vertex_attributes']) or \
    409                         (loaded_dict['vertex_attributes'] == None and \
     410        self.failUnless(num.alltrue(loaded_dict['vertex_attributes'] ==
     411                                    dict['vertex_attributes']) or
     412                        (loaded_dict['vertex_attributes'] == None and
    410413                         dict['vertex_attributes'] == []),
    411414                        fail_string + ' failed!! Test vertex_attributes')
    … …  
    416419        for seg, ldseg in map(None, dict['segment_tags'],
    417420                              loaded_dict['segment_tags']):
    418             self.failUnless(seg == ldseg,
    419                             fail_string + ' failed!! Test 8')
     421            msg = ('seg=\n"%s"\nshould be same as ldseg=\n"%s"'
     422                   % (str(seg), str(ldseg)))
     423            self.failUnless(seg == ldseg, fail_string + ' failed\n' + msg)
    420424
    421425        dict_array = num.array(dict['vertex_attribute_titles'])
    … …  
    429433
    430434        try:
     435            msg = ("loaded_dict['geo_reference']=\n%s\n"
     436                   "should be same as dict['geo_reference']=%s"
     437                  % (str(loaded_dict['geo_reference']),
     438                     str(dict['geo_reference'])))
    431439            self.failUnless(loaded_dict['geo_reference'] ==
    432440                            dict['geo_reference'],
    433                             fail_string + ' failed!! Test geo_reference')
     441                            fail_string + ' failed\n' + msg)
    434442        except KeyError:        #??# 2 lines below??
    435 #           self.failUnless(not dict.has_key('geo_reference' and
    436 #                               loaded_dict['geo_reference'] == None),
    437 #                           fail_string + ' failed!! Test geo_reference')
     443            msg = ("'dict' has no key 'geo_reference' "
     444                   "but loaded_dict['geo_reference'] isn't None")
    438445            self.failUnless(not dict.has_key('geo_reference') and
    439446                                loaded_dict['geo_reference'] == None,
    440                             fail_string + ' failed!! Test geo_reference')
     447                            fail_string + ' failed\n' + msg)
    441448
    442449########################## BAD .MSH ##########################

  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6410 r6428  
    27952795        # Invoke interpolation for vertex points       
    27962796        points = num.concatenate( (x[:,num.newaxis],y[:,num.newaxis]), axis=1 )
     2797        points = num.ascontiguousarray(points)
    27972798        sww2pts(self.domain.get_name(),
    27982799                quantity = 'elevation',
    … …  
    1090810909        domain.set_quantity('xmomentum', uh)
    1090910910        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
    10910 
    1091110911        for t in domain.evolve(yieldstep=1, finaltime = t_end):
    1091210912            pass
    … …  
    1133711337if __name__ == "__main__":
    1133811338    suite = unittest.makeSuite(Test_Data_Manager,'test')
    11339     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_maximum_inundation')
    1134011339
    1134111340    # FIXME (Ole): This is the test that fails under Windows

  • branches/numpy/anuga/shallow_water/test_smf.py

    r6410 r6428  
    5252        slide = slide_tsunami(length=len, depth=dep, slope=th, x0=x0, \
    5353                              width = wid, thickness=thk, kappa=kappa, kappad=kappad, \
    54                               verbose=False)
     54                              verbose=False)
    5555
    5656        assert num.allclose(slide.a3D, 0.07775819)
    … …  
    111111        slide = slide_tsunami(length, dep, th, x0, y0, \
    112112                              wid, thk, kappa, kappad, \
    113                               domain=domain,verbose=False)
     113                              domain=domain,verbose=False)
    114114
    115115        domain.set_quantity('stage', slide)
    116         stage = domain.get_quantity('stage')
    117         w = stage.get_values()
     116        stage = domain.get_quantity('stage')
     117        w = stage.get_values()
    118118
    119 ##      check = [[-0.0 -0.0 -0.0],
     119##        check = [[-0.0 -0.0 -0.0],
    120120##                 [-.189709745 -517.877716 -0.0],
    121121##                 [-0.0 -0.0 -2.7695931e-08],
    … …  
    126126##                 [-0.0 -0.0 -0.0]]
    127127
    128         assert num.allclose(min(min(w)), -517.877771593)
    129         assert num.allclose(max(max(w)), 0.0)
     128        assert num.allclose(num.min(w), -517.877771593)
     129        assert num.allclose(num.max(w), 0.0)
    130130        assert num.allclose(slide.a3D, 518.38797486)
    131131

  • branches/numpy/anuga/utilities/numerical_tools.py

    r6304 r6428  
    284284
    285285    n = num.searchsorted(num.sort(a), bins)
    286     n = num.concatenate( [n, [len(a)]] )
     286    n = num.concatenate([n, [len(a)]], axis=0)    #??default#
    287287
    288288    hist = n[1:]-n[:-1]

  • branches/numpy/anuga/utilities/system_tools.py

    r6415 r6428  
    294294    '''Convert 1-D list of strings to 2-D list of chars.'''
    295295
     296    if not l:
     297        return []
     298
     299    if l == ['']:
     300        l = [' ']
     301
    296302    maxlen = reduce(max, map(len, l))
    297303    ll = [x.ljust(maxlen) for x in l]

  • branches/numpy/anuga/utilities/test_system_tools.py

    r6415 r6428  
    55import numpy as num
    66import zlib
     7import os
    78from os.path import join, split, sep
     9from Scientific.IO.NetCDF import NetCDFFile
    810from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    9 from anuga.config import netcdf_float
     11from anuga.config import netcdf_float, netcdf_char, netcdf_int
    1012
    1113
    … …  
    5759
    5860        os.remove(tmp_name)
    59        
     61
    6062
    6163
    … …  
    7476
    7577        os.remove(tmp_name)       
    76        
     78
    7779        # Binary NetCDF File X 2 (use mktemp's name)
    7880
    … …  
    9496            fid.variables['test_array'][:] = test_array
    9597            fid.close()
    96            
     98
    9799            # Second file
    98100            filename2 = mktemp(suffix='.nc', dir='.')
    … …  
    103105            fid.variables['test_array'][:] = test_array
    104106            fid.close()
    105            
    106            
     107
     108
    107109            checksum1 = compute_checksum(filename1)
    108110            checksum2 = compute_checksum(filename2)       
    … …  
    125127        if path == '':
    126128            path = '.' + sep
    127        
     129
    128130        filename = path + sep +  'crc_test_file.png'
    129131
    … …  
    190192
    191193        MAX_CHARS = 10
    192         MAX_ENTRIES = 100000
     194        MAX_ENTRIES = 10000
    193195        A_INT = ord('a')
    194196        Z_INT = ord('z')
    … …  
    208210        self.failUnlessEqual(new_str_list, str_list)
    209211
     212    # special test - input list is ['']
     213    def test_string_to_char2(self):
     214        # generate a special list shown bad in load_mesh testing
     215        str_list = ['']
     216
     217        x = string_to_char(str_list)
     218        new_str_list = char_to_string(x)
     219
     220        self.failUnlessEqual(new_str_list, str_list)
     221
    210222
    211223################################################################################
    … …  
    214226################################################################################
    215227
     228    ##
     229    # @brief Helper function to write a list of strings to a NetCDF file.
     230    # @param filename Path to the file to write.
     231    # @param l The list of strings to write.
     232    def helper_write_msh_file(self, filename, l):
     233        # open the NetCDF file
     234        fd = NetCDFFile(filename, netcdf_mode_w)
     235        fd.description = 'Test file - string arrays'
     236
     237        # convert list of strings to num.array
     238        al = num.array(string_to_char(l), num.character)
     239
     240        # write the list
     241        fd.createDimension('num_of_strings', al.shape[0])
     242        fd.createDimension('size_of_strings', al.shape[1])
     243
     244        var = fd.createVariable('strings', netcdf_char,
     245                                ('num_of_strings', 'size_of_strings'))
     246        var[:] = al
     247
     248        fd.close()
     249
     250
     251    ##
     252    # @brief Helper function to read a NetCDF file and return a list of strings.
     253    # @param filename Path to the file to read.
     254    # @return A list of strings from the file.
     255    def helper_read_msh_file(self, filename):
     256        fid = NetCDFFile(filename, netcdf_mode_r)
     257        mesh = {}
     258
     259        # Get the 'strings' variable
     260        strings = fid.variables['strings'][:]
     261
     262        fid.close()
     263
     264        return char_to_string(strings)
     265
     266
     267    # test random strings to a NetCDF file
     268    def test_string_to_netcdf(self):
     269        import random
     270
     271        MAX_CHARS = 10
     272        MAX_ENTRIES = 10000
     273
     274        A_INT = ord('a')
     275        Z_INT = ord('z')
     276
     277        FILENAME = 'test.msh'
     278
     279        # generate some random strings in a list, with guaranteed lengths
     280        str_list = ['x' * MAX_CHARS]        # make first maximum length
     281        for entry in xrange(MAX_ENTRIES):
     282            length = random.randint(1, MAX_CHARS)
     283            s = ''
     284            for c in range(length):
     285                s += chr(random.randint(A_INT, Z_INT))
     286            str_list.append(s)
     287
     288        self.helper_write_msh_file(FILENAME, str_list)
     289        new_str_list = self.helper_read_msh_file(FILENAME)
     290
     291        self.failUnlessEqual(new_str_list, str_list)
     292        os.remove(FILENAME)
     293
     294    # special test - list [''] to a NetCDF file
     295    def test_string_to_netcdf2(self):
     296        FILENAME = 'test.msh'
     297
     298        # generate some random strings in a list, with guaranteed lengths
     299        str_list = ['']
     300
     301        self.helper_write_msh_file(FILENAME, str_list)
     302        new_str_list = self.helper_read_msh_file(FILENAME)
     303
     304        self.failUnlessEqual(new_str_list, str_list)
     305        os.remove(FILENAME)
     306
     307#        print ('test_string_to_char: str_list=%s, new_str_list=%s'
     308                #               % (str(str_list), str(new_str_list)))
    216309################################################################################
    217310
    218311if __name__ == "__main__":
    219     #suite = unittest.makeSuite(Test_system_tools, 'test')
    220     suite = unittest.makeSuite(Test_system_tools, 'test_string_to_char')
     312    suite = unittest.makeSuite(Test_system_tools, 'test')
    221313    runner = unittest.TextTestRunner()
    222314    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.