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

numpy changes

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.