Changeset 3954


Ignore:
Timestamp:
Nov 9, 2006, 10:41:21 AM (16 years ago)
Author:
ole
Message:

More cleanup of triangle and node formats + better comments

Location:
anuga_core/source
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

    r3945 r3954  
    240240
    241241    def get_nodes(self, absolute=False):
    242         """Return all node coordinates ordered in an Nx2 array.
     242        """Return all nodes in mesh.
     243
     244        The nodes are ordered in an Nx2 array where N is the number of nodes.
    243245        This is the same format they were provided in the constructor
    244246        i.e. without any duplication.
     
    260262       
    261263
    262     def get_vertex_coordinates(self, unique=False, absolute=False):
    263         """Return all vertex coordinates.
    264         Return all vertex coordinates for all triangles as a 3*N x 2 array
    265         where the jth vertex of the ith triangle is located in row 3*i+j.
    266 
    267         Boolean keyword unique will cause the points to be returned as
    268         they were provided in the constructor i.e. without any duplication
    269         in an N x 2 array.
    270 
    271         """
    272 
    273         if unique is True:       
    274             return self.get_nodes(absolute)
    275 
    276            
    277         V = self.vertex_coordinates
     264    def get_vertex_coordinates(self, absolute=False):
     265        """Return vertex coordinates for all triangles.
     266       
     267        Return all vertex coordinates for all triangles as a 3*M x 2 array
     268        where the jth vertex of the ith triangle is located in row 3*i+j and
     269        M the number of triangles in the mesh.
     270
     271        Boolean keyword argument absolute determines whether coordinates
     272        are to be made absolute by taking georeference into account
     273        Default is False as many parts of ANUGA expects relative coordinates.
     274        """
     275
     276
     277       
     278        V = self.vertex_coordinates #[:3*M,:]
    278279        if absolute is True:
    279280            if not self.geo_reference.is_absolute():
     
    294295
    295296    def compute_vertex_coordinates(self):
    296         """Return all vertex coordinates for all triangles as a 3*N x 2 array
     297        """Return all vertex coordinates for all triangles as a 3*M x 2 array
    297298        where the jth vertex of the ith triangle is located in row 3*i+j.
    298299
     
    301302        """
    302303
    303         N = self.number_of_triangles
    304         vertex_coordinates = zeros((3*N, 2), Float)
    305 
    306         for i in range(N):
     304        M = self.number_of_triangles
     305        vertex_coordinates = zeros((3*M, 2), Float)
     306
     307        for i in range(M):
    307308            for j in range(3):
    308                 k = self.triangles[i,j]  #Index of vertex 0
    309                 v_k = self.nodes[k]
    310 
    311                 vertex_coordinates[3*i+j,:] = v_k
    312 
     309                k = self.triangles[i,j] # Index of vertex j in triangle i
     310                vertex_coordinates[3*i+j,:] = self.nodes[k]
    313311
    314312        return vertex_coordinates
    315313
    316314
    317     def get_vertices(self, indices=None):
    318         """Get connectivity
    319         indices is the set of element ids of interest
    320         """
    321 
    322         N = self.number_of_full_triangles
     315
     316    def get_triangles(self, indices=None):
     317        """Get mesh triangles.
     318
     319        Return Mx3 integer array where M is the number of triangles.
     320        Each row corresponds to one triangle and the three entries are
     321        indices into the mesh nodes which can be obtained using the method
     322        get_nodes()
     323
     324        Optional argument, indices is the set of triangle ids of interest.
     325        """
     326
     327        M = self.number_of_full_triangles
    323328
    324329        if indices is None:
    325             #indices = range(len(self))  #len(self)=number of elements
    326             indices = range(N)
    327 
    328         return  take(self.triangles, indices)
     330            indices = range(M)
     331
     332        return take(self.triangles, indices)
    329333   
    330334
    331     #FIXME - merge these two (get_vertices and get_triangles)
    332     def get_triangles(self, obj=False):
    333         """Get connetivity
    334         Return triangles (triplets of indices into point coordinates)
    335        
    336         If obj is True return structure commensurate with replicated
    337         points, allowing for discontinuities
    338         (FIXME: Need good name for this concept)
    339         """
    340 
    341         if obj is True:
    342             m = len(self)  #Number of triangles
    343             M = 3*m        #Total number of unique vertices
    344             T = reshape(array(range(M)).astype(Int), (m,3))
    345         else:
    346             T = self.triangles
     335
     336    def get_disconnected_triangles(self):
     337        """Get mesh based on nodes obtained from get_vertex_coordinates.
     338
     339        Return array Mx3 array of integers where each row corresponds to
     340        a triangle. A triangle is a triplet of indices into
     341        point coordinates obtained from get_vertex_coordinates and each
     342        index appears only once
     343
     344        This provides a mesh where no triangles share nodes
     345        (hence the name disconnected triangles) and different
     346        nodes may have the same coordinates.
     347
     348        This version of the mesh is useful for storing meshes with
     349        discontinuities at each node and is e.g. used for storing
     350        data in sww files.
     351
     352        The triangles created will have the format
     353
     354        [[0,1,2],
     355         [3,4,5],
     356         [6,7,8],
     357         ...
     358         [3*M-3 3*M-2 3*M-1]]         
     359        """
     360
     361        M = len(self) # Number of triangles
     362        K = 3*M       # Total number of unique vertices
     363        T = reshape(array(range(K)).astype(Int), (M,3))
    347364
    348365        return T     
     
    351368
    352369    def get_unique_vertices(self,  indices=None):
    353         triangles = self.get_vertices(indices=indices)
     370        """FIXME(Ole): This function needs a docstring
     371        """
     372        triangles = self.get_triangles(indices=indices)
    354373        unique_verts = {}
    355374        for triangle in triangles:
     
    361380
    362381    def build_vertexlist(self):
    363         """Build vertexlist index by vertex ids and for each entry (point id)
     382        """Build vertexlist indexed by vertex ids and for each entry (point id)
    364383        build a list of (triangles, vertex_id) pairs that use the point
    365384        as vertex.
    366385
     386        The vertex list will have length N, where N is the number of nodes
     387        in the mesh.
     388
    367389        Preconditions:
    368390          self.nodes and self.triangles are defined
     
    372394        """
    373395
    374         vertexlist = [None]*len(self.nodes)
     396        vertexlist = [None]*self.number_of_nodes
    375397        for i in range(self.number_of_triangles):
    376398
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r3945 r3954  
    1515"""
    1616
    17 from Numeric import array, zeros, Float, less, concatenate, NewAxis, argmax, allclose
     17from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
     18     argmax, allclose
    1819from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
    1920
     
    10261027
    10271028
    1028     #Method for outputting model results
    1029     #FIXME: Split up into geometric and numeric stuff.
    1030     #FIXME: Geometric (X,Y,V) should live in mesh.py
    1031     #FIXME: STill remember to move XY to mesh
     1029    # Methods for outputting model results
    10321030    def get_vertex_values(self,
    10331031                          xy=True,
    1034                           smooth = None,
    1035                           precision = None,
    1036                           reduction = None):
    1037         """Return vertex values like an OBJ format
     1032                          smooth=None,
     1033                          precision=None,
     1034                          reduction=None):
     1035        """Return vertex values like an OBJ format i.e. one value per node.
    10381036
    10391037        The vertex values are returned as one sequence in the 1D float array A.
     
    10411039
    10421040        The connectivity is represented as an integer array, V, of dimension
    1043         M x 3, where M is the number of volumes. Each row has three indices
    1044         into the X, Y, A arrays defining the triangle.
     1041        Mx3, where M is the number of triangles. Each row has three indices
     1042        defining the triangle and they correspond to elements in the arrays
     1043        X, Y and A.
    10451044
    10461045        if smooth is True, vertex values corresponding to one common
    10471046        coordinate set will be smoothed according to the given
    10481047        reduction operator. In this case vertex coordinates will be
    1049         de-duplicated.
     1048        de-duplicated corresponding to the original nodes as obtained from
     1049        the method general_mesh.get_nodes()
    10501050
    10511051        If no smoothings is required, vertex coordinates and values will
    10521052        be aggregated as a concatenation of values at
    1053         vertices 0, vertices 1 and vertices 2
     1053        vertices 0, vertices 1 and vertices 2. This corresponds to
     1054        the node coordinates obtained from the method
     1055        general_mesh.get_vertex_coordinates()
    10541056
    10551057
     
    10661068
    10671069        if smooth is None:
     1070            # Take default from domain
    10681071            smooth = self.domain.smooth
    10691072
    10701073        if precision is None:
    10711074            precision = Float
    1072 
    1073         #Create connectivity
    1074 
    1075         if smooth == True:
     1075           
     1076
     1077        if smooth is True:
     1078            # Ensure continuous vertex values by combining
     1079            # values at each node using the reduction operator
    10761080           
    10771081            if reduction is None:
     1082                # Take default from domain               
    10781083                reduction = self.domain.reduction
    10791084
    1080             V = self.domain.get_vertices()
    1081             N = len(self.domain.vertexlist)
     1085            V = self.domain.get_triangles()
     1086            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    10821087            A = zeros(N, precision)
    1083 
    1084             #Smoothing loop
     1088            points = self.domain.get_nodes()           
     1089
     1090            # Reduction loop
    10851091            for k in range(N):
    10861092                L = self.domain.vertexlist[k]
    10871093
    1088                 #Go through all triangle, vertex pairs
    1089                 #contributing to vertex k and register vertex value
    1090 
    1091                 if L is None: continue #In case there are unused points
    1092 
     1094                # Go through all triangle, vertex pairs
     1095                # contributing to vertex k and register vertex value
     1096                if L is None: continue # In case there are unused points
    10931097                contributions = []
    10941098                for volume_id, vertex_id in L:
     
    10981102                A[k] = reduction(contributions)
    10991103
    1100 
    1101             if xy is True:
    1102                 X = self.domain.get_nodes()[:,0].astype(precision)
    1103                 Y = self.domain.get_nodes()[:,1].astype(precision)
    1104 
    1105                 return X, Y, A, V
    1106             else:
    1107                 return A, V
    1108         else:
    1109             #Don't smooth
    1110             #obj machinery moved to general_mesh
    1111 
    1112             # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
    1113             # These vert_id's will relate to the verts created below
    1114             #m = len(self.domain)  #Number of volumes
    1115             #M = 3*m        #Total number of unique vertices
    1116             #V = reshape(array(range(M)).astype(Int), (m,3))
    1117 
    1118             V = self.domain.get_triangles(obj=True)
    1119             #FIXME use get_vertices, when ready
    1120 
    1121             A = self.vertex_values.flat
    1122 
    1123             #Do vertex coordinates
    1124             if xy is True:
    1125                 C = self.domain.get_vertex_coordinates()
    1126 
    1127                 X = C[:,0:6:2].copy()
    1128                 Y = C[:,1:6:2].copy()
    1129 
    1130                 return X.flat, Y.flat, A, V
    1131             else:
    1132                 return A, V
     1104        else:
     1105            # Allow discontinuous vertex values
     1106            V = self.domain.get_disconnected_triangles()
     1107            points = self.domain.get_vertex_coordinates()
     1108            A = self.vertex_values.flat.astype(precision)
     1109
     1110
     1111        # Return   
     1112        if xy is True:
     1113            X = points[:,0].astype(precision)
     1114            Y = points[:,1].astype(precision)
     1115           
     1116            return X, Y, A, V
     1117        else:
     1118            return A, V           
     1119
    11331120
    11341121
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r3945 r3954  
    2525
    2626        #Create basic mesh
    27         points, vertices, boundary = rectangular(1, 3)
    28         domain = General_mesh(points, vertices)
     27        nodes, triangles, _ = rectangular(1, 3)
     28        domain = General_mesh(nodes, triangles)
    2929
    3030
    31         assert allclose(domain.get_vertex_coordinates(unique=True),
    32                         domain.nodes)
     31        assert allclose(domain.get_nodes(), nodes)
    3332
    34         #assert allclose(domain.get_vertex_coordinates(), ...TODO
    35         #assert allclose(domain.get_vertex_coordinates(absolute=True), ...TODO
     33
     34        M = domain.number_of_triangles       
     35
     36        V = domain.get_vertex_coordinates()
     37        assert V.shape[0] == 3*M
     38
     39        for i in range(M):
     40            for j in range(3):
     41                k = triangles[i,j]  #Index of vertex j in triangle i
     42                assert allclose(V[3*i+j,:], nodes[k])
     43
     44
    3645       
    3746       
     
    4453
    4554        #Create basic mesh
    46         points, vertices, boundary = rectangular(1, 3)
    47         domain = General_mesh(points, vertices)
     55        nodes, triangles, _ = rectangular(1, 3)
     56        domain = General_mesh(nodes, triangles)
    4857
    4958        value = [7]
    50         #indexes = [1]  #FIXME (Ole): Should this be used
    51         assert  domain.get_vertices() == domain.triangles
    52         assert domain.get_vertices([0,4]) == [domain.triangles[0],
    53                                               domain.triangles[4]]
     59        assert allclose(domain.get_triangles(), triangles)
     60        assert allclose(domain.get_triangles([0,4]),
     61                        [triangles[0], triangles[4]])
    5462       
    5563    def test_areas(self):
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r3953 r3954  
    372372                                      precision=self.precision)
    373373
    374 
    375 
    376 
    377         x[:] = X.astype(self.precision)
    378        
    379         y[:] = Y.astype(self.precision)
    380 
    381         # FIXME (HACK)
    382         if len(z) <> len(Z):
    383             truncation = self.domain.number_of_full_nodes       
    384             Z = Z[:truncation]
    385             print len(z), len(Z), truncation, len(self.domain.get_nodes())
    386        
    387         z[:] = Z.astype(self.precision)
     374        volumes[:] = V.astype(volumes.typecode())
     375        x[:] = X
     376        y[:] = Y
     377        z[:] = Z
    388378
    389379        #FIXME: Backwards compatibility
    390380        z = fid.variables['z']
    391         z[:] = Z.astype(self.precision)
     381        z[:] = Z
    392382        ################################
    393383
    394         volumes[:] = V.astype(volumes.typecode())
     384
    395385
    396386        #Close
     
    501491                                          precision = self.precision)
    502492
    503                 # HACK
    504                 truncation = self.domain.number_of_full_nodes
    505                 # HACK
    506                 if stage.shape[1] <> len(A):               
    507                     A = A[:truncation]
    508                
    509493
    510494                #FIXME: Make this general (see below)
     
    29762960    outfile.variables['x'][:] = x - xllcorner
    29772961    outfile.variables['y'][:] = y - yllcorner
    2978     outfile.variables['z'][:] = z             #FIXME HACK
     2962    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
    29792963    outfile.variables['elevation'][:] = z
    29802964    outfile.variables['time'][:] = times   #Store time relative
    2981     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
     2965    outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
    29822966
    29832967
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r3928 r3954  
    191191        value per vertex using self.reduction.
    192192        """
     193
     194        # FIXME (Ole): how about using the word continuous vertex values?
    193195        self.smooth = not flag
    194196
  • anuga_core/source/anuga_parallel/print_stats.py

    r3514 r3954  
    4646def build_full_flag(domain, ghost_recv_dict):
    4747
    48     tri_full_flag = ones(len(domain.get_vertices()), Int8)
     48    tri_full_flag = ones(len(domain.get_triangles()), Int8)
    4949    for i in ghost_recv_dict.keys():
    5050        for id in ghost_recv_dict[i][0]:
Note: See TracChangeset for help on using the changeset viewer.