Changeset 1237


Ignore:
Timestamp:
Apr 20, 2005, 5:41:03 PM (20 years ago)
Author:
steve
Message:
 
Location:
inundation/ga/storm_surge
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution-parallel/domain.py

    r1232 r1237  
    430430        self.distribute_to_vertices_and_edges()
    431431
     432        #Initialize boundary values
     433        self.update_boundary()
    432434
    433435        #Or maybe restore from latest checkpoint
    … …  
    439441
    440442        while True:
     443
     444            #Compute fluxes across each element edge
     445            self.compute_fluxes()
     446
     447            #Update timestep to fit yieldstep and finaltime
     448            self.update_timestep(yieldstep, finaltime)
     449
     450            #Update conserved quantities
     451            self.update_conserved_quantities()
     452
     453            #Update vertex and edge values
     454            self.distribute_to_vertices_and_edges()
     455
    441456            #Update boundary values
    442457            self.update_boundary()
    443 
    444             #Compute fluxes across each element edge
    445             self.compute_fluxes()
    446 
    447             #Update timestep to fit yieldstep and finaltime
    448             self.update_timestep(yieldstep, finaltime)
    449 
    450             #Update conserved quantities
    451             self.update_conserved_quantities()
    452 
    453             #Update vertex and edge values
    454             self.distribute_to_vertices_and_edges()
    455458
    456459            #Update time

  • inundation/ga/storm_surge/pyvolution-parallel/general_mesh.py

    r1178 r1237  
    66
    77    A triangular element is defined in terms of three vertex ids,
    8     ordered counter clock-wise, 
     8    ordered counter clock-wise,
    99    each corresponding to a given coordinate set.
    1010    Vertices from different elements can point to the same
    … …  
    1212
    1313    Coordinate sets are implemented as an N x 2 Numeric array containing
    14     x and y coordinates. 
    15    
     14    x and y coordinates.
     15
    1616
    1717    To instantiate:
    … …  
    3333        c = [2.0,0.0]
    3434        e = [2.0, 2.0]
    35        
     35
    3636        points = [a, b, c, e]
    37         triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
     37        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
    3838        mesh = Mesh(points, triangles)
    39    
     39
    4040        #creates two triangles: bac and bce
    4141
    4242    Other:
    43    
     43
    4444      In addition mesh computes an Nx6 array called vertex_coordinates.
    45       This structure is derived from coordinates and contains for each 
     45      This structure is derived from coordinates and contains for each
    4646      triangle the three x,y coordinates at the vertices.
    47                                
     47
    4848
    4949        This is a cut down version of mesh from pyvolution\mesh.py
    … …  
    5858
    5959        origin is a 3-tuple consisting of UTM zone, easting and northing.
    60         If specified coordinates are assumed to be relative to this origin. 
     60        If specified coordinates are assumed to be relative to this origin.
    6161        """
    6262
    … …  
    6767
    6868        if geo_reference is None:
    69             self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) 
     69            self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0)
    7070        else:
    7171            self.geo_reference = geo_reference
    … …  
    7676
    7777        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
    78         assert len(self.coordinates.shape) == 2, msg       
     78        assert len(self.coordinates.shape) == 2, msg
    7979
    8080        msg = 'Vertex indices reference non-existing coordinate sets'
    … …  
    8484        #Register number of elements (N)
    8585        self.number_of_elements = N = self.triangles.shape[0]
    86    
    87         #Allocate space for geometric quantities       
     86
     87        #Allocate space for geometric quantities
    8888        self.normals = zeros((N, 6), Float)
    8989        self.areas = zeros(N, Float)
    … …  
    9292        #Get x,y coordinates for all triangles and store
    9393        self.vertex_coordinates = V = self.compute_vertex_coordinates()
    94        
     94
    9595        #Initialise each triangle
    9696        for i in range(N):
    9797            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
    98            
     98
    9999            x0 = V[i, 0]; y0 = V[i, 1]
    100100            x1 = V[i, 2]; y1 = V[i, 3]
    101             x2 = V[i, 4]; y2 = V[i, 5]           
     101            x2 = V[i, 4]; y2 = V[i, 5]
    102102
    103103            #Area
    … …  
    107107            msg += ' is degenerate:  area == %f' %self.areas[i]
    108108            assert self.areas[i] > 0.0, msg
    109        
     109
    110110
    111111            #Normals
    … …  
    120120
    121121            n0 = array([x2 - x1, y2 - y1])
    122             l0 = sqrt(sum(n0**2))
     122            l0 = sqrt(sum(n0**2))
    123123
    124124            n1 = array([x0 - x2, y0 - y2])
    … …  
    137137                                  n1[1], -n1[0],
    138138                                  n2[1], -n2[0]]
    139            
     139
    140140            #Edgelengths
    141141            self.edgelengths[i, :] = [l0, l1, l2]
    142142
    143143        #Build vertex list
    144         self.build_vertexlist()       
    145            
     144        self.build_vertexlist()
     145
    146146    def __len__(self):
    147147        return self.number_of_elements
    … …  
    154154        """Return all normal vectors.
    155155        Return normal vectors for all triangles as an Nx6 array
    156         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
     156        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
    157157        """
    158158        return self.normals
    … …  
    163163        Return value is the numeric array slice [x, y]
    164164        """
    165         return self.normals[i, 2*j:2*j+2]     
    166    
    167 
    168            
     165        return self.normals[i, 2*j:2*j+2]
     166
     167
     168
    169169    def get_vertex_coordinates(self):
    170170        """Return all vertex coordinates.
    171171        Return all vertex coordinates for all triangles as an Nx6 array
    172         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
     172        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
    173173        """
    174174        return self.vertex_coordinates
    … …  
    179179        Return value is the numeric array slice [x, y]
    180180        """
    181         return self.vertex_coordinates[i, 2*j:2*j+2]     
     181        return self.vertex_coordinates[i, 2*j:2*j+2]
    182182
    183183
    … …  
    189189        #FIXME: Perhaps they should be ordered as in obj files??
    190190        #See quantity.get_vertex_values
    191        
     191
    192192        from Numeric import zeros, Float
    193193
    194194        N = self.number_of_elements
    195         vertex_coordinates = zeros((N, 6), Float) 
     195        vertex_coordinates = zeros((N, 6), Float)
    196196
    197197        for i in range(N):
    … …  
    203203
    204204        return vertex_coordinates
    205    
     205
    206206    def get_vertices(self,  indexes=None):
    207207        """Get connectivity
    208208        indexes is the set of element ids of interest
    209209        """
    210        
     210
    211211        from Numeric import take
    212        
     212
    213213        if (indexes ==  None):
    214214            indexes = range(len(self))  #len(self)=number of elements
    215            
     215
    216216        return  take(self.triangles, indexes)
    217217
    … …  
    220220        unique_verts = {}
    221221        for triangle in triangles:
    222            unique_verts[triangle[0]] = 0 
    223            unique_verts[triangle[1]] = 0 
     222           unique_verts[triangle[0]] = 0
     223           unique_verts[triangle[1]] = 0
    224224           unique_verts[triangle[2]] = 0
    225         return unique_verts.keys()   
    226        
     225        return unique_verts.keys()
     226
    227227    def build_vertexlist(self):
    228228        """Build vertexlist index by vertex ids and for each entry (point id)
    … …  
    231231
    232232        Preconditions:
    233           self.coordinates and self.triangles are defined 
    234        
     233          self.coordinates and self.triangles are defined
     234
    235235        Postcondition:
    236236          self.vertexlist is built
    … …  
    242242            a = self.triangles[i, 0]
    243243            b = self.triangles[i, 1]
    244             c = self.triangles[i, 2]           
     244            c = self.triangles[i, 2]
    245245
    246246            #Register the vertices v as lists of
    … …  
    251251                    vertexlist[v] = []
    252252
    253                 vertexlist[v].append( (i, vertex_id) )   
     253                vertexlist[v].append( (i, vertex_id) )
    254254
    255255        self.vertexlist = vertexlist
    256    
     256
    257257
    258258    def get_extent(self):
    … …  
    264264        X = C[:,0:6:2].copy()
    265265        Y = C[:,1:6:2].copy()
    266        
     266
    267267        xmin = min(X.flat)
    268268        xmax = max(X.flat)
    269269        ymin = min(Y.flat)
    270         ymax = max(Y.flat)               
    271        
     270        ymax = max(Y.flat)
     271
    272272        return xmin, xmax, ymin, ymax
    273273
    … …  
    275275    def get_area(self):
    276276        """Return total area of mesh
    277         """
    278 
    279         return sum(self.areas)         
    280        
     277        """
     278
     279        return sum(self.areas)
     280

  • inundation/ga/storm_surge/pyvolution-parallel/mesh.py

    r1193 r1237  
    6767        """
    6868
    69        
     69
    7070
    7171        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
    … …  
    113113                #a triangle to the midpoint of the side of the triangle
    114114                #closest to the centroid
    115                 d0 = sqrt(sum( (centroid-m0)**2 ))
     115                d0 = sqrt(sum( (centroid-m0)**2 ))
    116116                d1 = sqrt(sum( (centroid-m1)**2 ))
    117117                d2 = sqrt(sum( (centroid-m2)**2 ))
    … …  
    129129
    130130                self.radii[i]=self.areas[i]/(2*(a+b+c))
    131                
     131
    132132
    133133            #Initialise Neighbours (-1 means that it is a boundary neighbour)
    … …  
    446446            #Note: Could be arbitrary, but nice to have
    447447            #a unique way of selecting
    448             dist_A = sqrt(sum( (A-pmin)**2 ))
    449             dist_B = sqrt(sum( (B-pmin)**2 ))
     448            dist_A = sqrt(sum( (A-pmin)**2 ))
     449            dist_B = sqrt(sum( (B-pmin)**2 ))
    450450
    451451            #Find minimal point
    … …  
    584584        #    assert self.neighbours[id,edge] < 0
    585585        #
    586         #NOTE (Ole): I reckon this was resolved late 2004?
    587         #
    588         #See domain.set_boundary
     586        #NOTE (Ole): I reckon this was resolved late 2004?
     587        #
     588        #See domain.set_boundary
    589589
    590590

  • inundation/ga/storm_surge/pyvolution-parallel/quantity.py

    r1011 r1237  
    66
    77   domain: Associated domain structure. Required.
    8    
     8
    99   vertex_values: N x 3 array of values at each vertex for each element.
    1010                  Default None
    … …  
    2828
    2929        if vertex_values is None:
    30             N = domain.number_of_elements           
     30            N = domain.number_of_elements
    3131            self.vertex_values = zeros((N, 3), Float)
    32         else:   
     32        else:
    3333            self.vertex_values = array(vertex_values).astype(Float)
    3434
    … …  
    3636            assert V == 3,\
    3737                   'Three vertex values per element must be specified'
    38            
     38
    3939
    4040            msg = 'Number of vertex values (%d) must be consistent with'\
    … …  
    4242            msg += 'number of elements in specified domain (%d).'\
    4343                   %domain.number_of_elements
    44            
     44
    4545            assert N == domain.number_of_elements, msg
    4646
    47         self.domain = domain   
     47        self.domain = domain
    4848
    4949        #Allocate space for other quantities
    5050        self.centroid_values = zeros(N, Float)
    51         self.edge_values = zeros((N, 3), Float)       
     51        self.edge_values = zeros((N, 3), Float)
    5252
    5353        #Intialise centroid and edge_values
    … …  
    5656    def __len__(self):
    5757        return self.centroid_values.shape[0]
    58    
     58
    5959    def interpolate(self):
    6060        """Compute interpolated values at edges and centroid
    … …  
    6262        """
    6363
    64         N = self.vertex_values.shape[0]       
     64        N = self.vertex_values.shape[0]
    6565        for i in range(N):
    6666            v0 = self.vertex_values[i, 0]
    6767            v1 = self.vertex_values[i, 1]
    6868            v2 = self.vertex_values[i, 2]
    69            
     69
    7070            self.centroid_values[i] = (v0 + v1 + v2)/3
    7171
    … …  
    7777        #(either from this module or C-extension)
    7878        interpolate_from_vertices_to_edges(self)
    79            
    80            
     79
     80
     81
    8182    def set_values(self, X, location='vertices', indexes = None):
    8283        """Set values for quantity
    … …  
    9697        If values are described a function, it will be evaluated at
    9798        specified points
    98        
     99
    99100        If indexex is not 'unique vertices' Indexes is the set of element ids
    100101        that the operation applies to.
    101102        If indexex is 'unique vertices' Indexes is the set of vertex ids
    102103        that the operation applies to.
    103        
    104        
     104
     105
    105106        If selected location is vertices, values for centroid and edges
    106107        will be assigned interpolated values.
    … …  
    115116        if X is None:
    116117            msg = 'Given values are None'
    117             raise msg           
     118            raise msg
    118119
    119120        import types, Numeric
    … …  
    122123                                 'Indices must be a list or None'
    123124
    124        
     125
    125126        if callable(X):
    126127            #Use function specific method
    127             self.set_function_values(X, location, indexes = indexes)       
     128            self.set_function_values(X, location, indexes = indexes)
    128129        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
    129130            if location == 'centroids':
    … …  
    134135                    for i in indexes:
    135136                        self.centroid_values[i,:] = X
    136                        
     137
    137138            elif location == 'edges':
    138139                if (indexes ==  None):
    … …  
    142143                    for i in indexes:
    143144                        self.edge_values[i,:] = X
    144                        
     145
    145146            elif location == 'unique vertices':
    146147                if (indexes ==  None):
    … …  
    151152                    for unique_vert_id in indexes:
    152153                        triangles = self.domain.vertexlist[unique_vert_id]
    153                        
     154
    154155                        #In case there are unused points
    155                         if triangles is None: continue 
     156                        if triangles is None: continue
    156157
    157158                        #Go through all triangle, vertex pairs
    … …  
    159160                        for triangle_id, vertex_id in triangles:
    160161                            self.vertex_values[triangle_id, vertex_id] = X
    161                
     162
    162163                        #Intialise centroid and edge_values
    163164                        self.interpolate()
    … …  
    168169                    #Brute force
    169170                    for i_vertex in indexes:
    170                         self.vertex_values[i_vertex,:] = X     
     171                        self.vertex_values[i_vertex,:] = X
    171172
    172173        else:
    … …  
    177178            #Intialise centroid and edge_values
    178179            self.interpolate()
    179            
     180
    180181        if location == 'centroids':
    181182            #Extrapolate 1st order - to capture notion of area being specified
    182             self.extrapolate_first_order()           
    183          
     183            self.extrapolate_first_order()
     184
     185
    184186    def get_values(self, location='vertices', indexes = None):
    185187        """get values for quantity
    … …  
    197199        (N if indexes = None).  Each value will be a list of the three
    198200        vertex values for this quantity.
    199        
     201
    200202        Indexes is the set of element ids that the operation applies to.
    201        
     203
    202204        """
    203205        from Numeric import take
    … …  
    205207        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
    206208            msg = 'Invalid location: %s' %location
    207             raise msg         
    208        
     209            raise msg
     210
    209211        import types, Numeric
    210212        assert type(indexes) in [types.ListType, types.NoneType,
    211213                                 Numeric.ArrayType],\
    212214                                 'Indices must be a list or None'
    213            
     215
    214216        if location == 'centroids':
    215217            if (indexes ==  None):
    216218                indexes = range(len(self))
    217             return take(self.centroid_values,indexes)                       
     219            return take(self.centroid_values,indexes)
    218220        elif location == 'edges':
    219221            if (indexes ==  None):
    220222                indexes = range(len(self))
    221             return take(self.edge_values,indexes)                     
     223            return take(self.edge_values,indexes)
    222224        elif location == 'unique vertices':
    223225            if (indexes ==  None):
    … …  
    227229            for unique_vert_id in indexes:
    228230                triangles = self.domain.vertexlist[unique_vert_id]
    229                        
     231
    230232                #In case there are unused points
    231                 if triangles is None: 
    232                     msg = 'Unique vertex not associated with triangles' 
    233                     raise msg 
     233                if triangles is None:
     234                    msg = 'Unique vertex not associated with triangles'
     235                    raise msg
    234236
    235237                # Go through all triangle, vertex pairs
    … …  
    280282            raise 'Not implemented: %s' %location
    281283
    282            
    283284    def set_array_values(self, values, location='vertices', indexes = None):
    284285        """Set values for quantity
    … …  
    288289        Permissible options are: vertices, edges, centroid, unique vertices
    289290        Default is "vertices"
    290        
     291
    291292        indexes - if this action is carried out on a subset of
    292293        elements or unique vertices
    293294        The element/unique vertex indexes are specified here.
    294        
     295
    295296        In case of location == 'centroid' the dimension values must
    296297        be a list of a Numerical array of length N, N being the number
    … …  
    315316            indexes = array(indexes).astype(Int)
    316317            msg = 'Number of values must match number of indexes'
    317             assert values.shape[0] == indexes.shape[0], msg 
    318            
     318            assert values.shape[0] == indexes.shape[0], msg
     319
    319320        N = self.centroid_values.shape[0]
    320        
     321
    321322        if location == 'centroids':
    322323            assert len(values.shape) == 1, 'Values array must be 1d'
    … …  
    325326                msg = 'Number of values must match number of elements'
    326327                assert values.shape[0] == N, msg
    327            
     328
    328329                self.centroid_values = values
    329330            else:
    330331                msg = 'Number of values must match number of indexes'
    331332                assert values.shape[0] == indexes.shape[0], msg
    332                
     333
    333334                #Brute force
    334335                for i in range(len(indexes)):
    335336                    self.centroid_values[indexes[i]] = values[i]
    336                    
     337
    337338        elif location == 'edges':
    338339            assert len(values.shape) == 2, 'Values array must be 2d'
    … …  
    340341            msg = 'Number of values must match number of elements'
    341342            assert values.shape[0] == N, msg
    342            
    343             msg = 'Array must be N x 3'           
     343
     344            msg = 'Array must be N x 3'
    344345            assert values.shape[1] == 3, msg
    345            
     346
    346347            self.edge_values = values
    347            
     348
    348349        elif location == 'unique vertices':
    349350            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
    … …  
    360361                 #   self.set_vertex_values(values)
    361362                #else:
    362                 #    for element_index, value in map(None, indexes, values): 
     363                #    for element_index, value in map(None, indexes, values):
    363364                #        self.vertex_values[element_index, :] = value
    364                        
     365
    365366            elif len(values.shape) == 2:
    366367                #Vertex values are given as a triplet for each triangle
    367                
    368                 msg = 'Array must be N x 3'           
     368
     369                msg = 'Array must be N x 3'
    369370                assert values.shape[1] == 3, msg
    370                
     371
    371372                if indexes == None:
    372373                    self.vertex_values = values
    373374                else:
    374                     for element_index, value in map(None, indexes, values): 
     375                    for element_index, value in map(None, indexes, values):
    375376                        self.vertex_values[element_index] = value
    376             else:   
     377            else:
    377378                msg = 'Values array must be 1d or 2d'
    378379                raise msg
    379            
    380 
    381     # FIXME have a get_vertex_values as well, so the 'stage' quantity can be
    382     # set, based on the elevation
     380
     381
     382        #FIXME have a get_vertex_values as well, so the 'stage' quantity can be
     383        # set, based on the elevation
     384
    383385    def set_vertex_values(self, A, indexes = None):
    384386        """Set vertex values for all unique vertices based on input array A
    … …  
    387389        one value for each row in vertexlist.
    388390
    389         indexes is the list of vertex_id's that will be set. 
     391        indexes is the list of vertex_id's that will be set.
    390392
    391393        Note: Functions not allowed
    … …  
    409411        for i_index,unique_vert_id in enumerate(vertex_list):
    410412            triangles = self.domain.vertexlist[unique_vert_id]
    411                
     413
    412414            if triangles is None: continue #In case there are unused points
    413415
    … …  
    416418            for triangle_id, vertex_id in triangles:
    417419                self.vertex_values[triangle_id, vertex_id] = A[i_index]
    418                
     420
    419421        #Intialise centroid and edge_values
    420422        self.interpolate()
    421  
     423
    422424    def smooth_vertex_values(self, value_array='field_values',
    423425                             precision = None):
    … …  
    432434        from Numeric import concatenate, zeros, Float, Int, array, reshape
    433435
    434        
     436
    435437        A,V = self.get_vertex_values(xy=False,
    436438                                     value_array=value_array,
    … …  
    450452        elif value_array == 'conserved_quantities':
    451453            Volume.interpolate_conserved_quantities()
    452        
    453    
    454     #Method for outputting model results
    455     #FIXME: Split up into geometric and numeric stuff.
    456     #FIXME: Geometric (X,Y,V) should live in mesh.py
    457     #FIXME: STill remember to move XY to mesh
     454
     455
    458456    def get_vertex_values(self,
    459                           xy=True, 
     457                          xy=True,
    460458                          smooth = None,
    461459                          precision = None,
    … …  
    465463        The vertex values are returned as one sequence in the 1D float array A.
    466464        If requested the coordinates will be returned in 1D arrays X and Y.
    467        
     465
    468466        The connectivity is represented as an integer array, V, of dimension
    469467        M x 3, where M is the number of volumes. Each row has three indices
    … …  
    474472        reduction operator. In this case vertex coordinates will be
    475473        de-duplicated.
    476        
     474
    477475        If no smoothings is required, vertex coordinates and values will
    478476        be aggregated as a concatenation of values at
    479477        vertices 0, vertices 1 and vertices 2
    480478
    481        
     479
    482480        Calling convention
    483481        if xy is True:
    484482           X,Y,A,V = get_vertex_values
    485483        else:
    486            A,V = get_vertex_values         
    487            
    488         """
     484           A,V = get_vertex_values
     485        """
     486        #Method for outputting model results
     487        #FIXME: Split up into geometric and numeric stuff.
     488        #FIXME: Geometric (X,Y,V) should live in mesh.py
     489        #FIXME: STill remember to move XY to mesh
    489490
    490491        from Numeric import concatenate, zeros, Float, Int, array, reshape
    491 
    492492
    493493        if smooth is None:
    … …  
    499499        if reduction is None:
    500500            reduction = self.domain.reduction
    501            
     501
    502502        #Create connectivity
    503                        
     503
    504504        if smooth == True:
    505505
    … …  
    511511            for k in range(N):
    512512                L = self.domain.vertexlist[k]
    513                
     513
    514514                #Go through all triangle, vertex pairs
    515515                #contributing to vertex k and register vertex value
    516516
    517517                if L is None: continue #In case there are unused points
    518                
     518
    519519                contributions = []
    520520                for volume_id, vertex_id in L:
    … …  
    522522                    contributions.append(v)
    523523
    524                 A[k] = reduction(contributions)   
     524                A[k] = reduction(contributions)
    525525
    526526
    … …  
    528528                X = self.domain.coordinates[:,0].astype(precision)
    529529                Y = self.domain.coordinates[:,1].astype(precision)
    530                
     530
    531531                return X, Y, A, V
    532532            else:
    … …  
    540540            M = 3*m        #Total number of unique vertices
    541541            V = reshape(array(range(M)).astype(Int), (m,3))
    542            
     542
    543543            A = self.vertex_values.flat
    544544
    545             #Do vertex coordinates   
    546             if xy is True:               
     545            #Do vertex coordinates
     546            if xy is True:
    547547                C = self.domain.get_vertex_coordinates()
    548548
    549549                X = C[:,0:6:2].copy()
    550                 Y = C[:,1:6:2].copy()               
    551 
    552                 return X.flat, Y.flat, A, V           
     550                Y = C[:,1:6:2].copy()
     551
     552                return X.flat, Y.flat, A, V
    553553            else:
    554554                return A, V
    555555
    556            
    557556    def extrapolate_first_order(self):
    558557        """Extrapolate conserved quantities from centroid to
    … …  
    560559        first order scheme.
    561560        """
    562        
     561
    563562        qc = self.centroid_values
    564563        qv = self.vertex_values
    … …  
    566565        for i in range(3):
    567566            qv[:,i] = qc
    568            
     567
    569568
    570569    def get_integral(self):
    571570        """Compute the integral of quantity across entire domain
    572         """
    573         integral = 0
    574         for k in range(self.domain.number_of_elements):
     571        """
     572        integral = 0
     573        for k in range(self.domain.number_of_elements):
    575574            area = self.domain.areas[k]
    576             qc = self.centroid_values[k]     
     575            qc = self.centroid_values[k]
    577576            integral += qc*area
    578            
    579         return integral   
    580        
     577
     578        return integral
     579
    581580
    582581class Conserved_quantity(Quantity):
    … …  
    590589    def __init__(self, domain, vertex_values=None):
    591590        Quantity.__init__(self, domain, vertex_values)
    592        
     591
    593592        from Numeric import zeros, Float
    594593
    … …  
    603602        self.explicit_update = zeros(N, Float )
    604603        self.semi_implicit_update = zeros(N, Float )
    605                
     604
    606605
    607606    def update(self, timestep):
    … …  
    615614        #(either from this module or C-extension)
    616615        return compute_gradients(self)
    617            
     616
    618617
    619618    def limit(self):
    … …  
    621620        #(either from this module or C-extension)
    622621        limit(self)
    623        
    624        
     622
     623
    625624    def extrapolate_second_order(self):
    626625        #Call correct module function
    … …  
    629628
    630629
    631 def update(quantity, timestep):   
     630def update(quantity, timestep):
    632631    """Update centroid values based on values stored in
    633632    explicit_update and semi_implicit_update as well as given timestep
    … …  
    642641        domain.quantities['ymomentum'].explicit_update = ...
    643642
    644        
     643
    645644
    646645    Explicit terms must have the form
    647    
     646
    648647        G(q, t)
    649    
     648
    650649    and explicit scheme is
    651    
     650
    652651       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
    653652
    654653
    655654    Semi implicit forcing terms are assumed to have the form
    656    
     655
    657656       G(q, t) = H(q, t) q
    658    
     657
    659658    and the semi implicit scheme will then be
    660    
     659
    661660      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
    662661
    663    
     662
    664663    """
    665    
     664
    666665    from Numeric import sum, equal, ones, Float
    667        
     666
    668667    N = quantity.centroid_values.shape[0]
    669668
    … …  
    673672
    674673    for k in range(N):
    675         x = quantity.centroid_values[k] 
     674        x = quantity.centroid_values[k]
    676675        if x == 0.0:
    677676            #FIXME: Is this right
    678             quantity.semi_implicit_update[k] = 0.0           
    679         else:
    680             quantity.semi_implicit_update[k] /= x             
    681            
     677            quantity.semi_implicit_update[k] = 0.0
     678        else:
     679            quantity.semi_implicit_update[k] /= x
     680
    682681    #Explicit updates
    683682    quantity.centroid_values += timestep*quantity.explicit_update
    684            
     683
    685684    #Semi implicit updates
    686685    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
    … …  
    702701        q1 = quantity.vertex_values[k, 1]
    703702        q2 = quantity.vertex_values[k, 2]
    704            
     703
    705704        quantity.edge_values[k, 0] = 0.5*(q1+q2)
    706         quantity.edge_values[k, 1] = 0.5*(q0+q2) 
     705        quantity.edge_values[k, 1] = 0.5*(q0+q2)
    707706        quantity.edge_values[k, 2] = 0.5*(q0+q1)
    708707
    709708
    710709
    711 def extrapolate_second_order(quantity):       
     710def extrapolate_second_order(quantity):
    712711    """Extrapolate conserved quantities from centroid to
    713712    vertices for each volume using
    714713    second order scheme.
    715714    """
    716        
     715
    717716    a, b = quantity.compute_gradients()
    718717
    … …  
    720719    qc = quantity.centroid_values
    721720    qv = quantity.vertex_values
    722    
     721
    723722    #Check each triangle
    724723    for k in range(quantity.domain.number_of_elements):
    725         #Centroid coordinates           
     724        #Centroid coordinates
    726725        x, y = quantity.domain.centroid_coordinates[k]
    727        
     726
    728727        #vertex coordinates
    729728        x0, y0, x1, y1, x2, y2 = X[k,:]
    730        
     729
    731730        #Extrapolate
    732731        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
    733732        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
    734         qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)           
    735 
    736 
    737 def compute_gradients(quantity):                   
     733        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
     734
     735
     736def compute_gradients(quantity):
    738737    """Compute gradients of triangle surfaces defined by centroids of
    739738    neighbouring volumes.
    … …  
    744743    from Numeric import zeros, Float
    745744    from util import gradient
    746    
     745
    747746    centroid_coordinates = quantity.domain.centroid_coordinates
    748     surrogate_neighbours = quantity.domain.surrogate_neighbours   
    749     centroid_values = quantity.centroid_values   
    750     number_of_boundaries = quantity.domain.number_of_boundaries     
    751    
     747    surrogate_neighbours = quantity.domain.surrogate_neighbours
     748    centroid_values = quantity.centroid_values
     749    number_of_boundaries = quantity.domain.number_of_boundaries
     750
    752751    N = centroid_values.shape[0]
    753752
    754753    a = zeros(N, Float)
    755754    b = zeros(N, Float)
    756    
     755
    757756    for k in range(N):
    758757        if number_of_boundaries[k] < 2:
    759758            #Two or three true neighbours
    760759
    761             #Get indices of neighbours (or self when used as surrogate)     
     760            #Get indices of neighbours (or self when used as surrogate)
    762761            k0, k1, k2 = surrogate_neighbours[k,:]
    763762
    764             #Get data       
     763            #Get data
    765764            q0 = centroid_values[k0]
    766765            q1 = centroid_values[k1]
    767             q2 = centroid_values[k2]                   
     766            q2 = centroid_values[k2]
    768767
    769768            x0, y0 = centroid_coordinates[k0] #V0 centroid
    … …  
    773772            #Gradient
    774773            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    775        
     774
    776775        elif number_of_boundaries[k] == 2:
    777776            #One true neighbour
    … …  
    789788
    790789            x0, y0 = centroid_coordinates[k0] #V0 centroid
    791             x1, y1 = centroid_coordinates[k1] #V1 centroid       
     790            x1, y1 = centroid_coordinates[k1] #V1 centroid
    792791
    793792            #Gradient
    … …  
    798797
    799798        else:
    800             #No true neighbours -       
     799            #No true neighbours -
    801800            #Fall back to first order scheme
    802801            pass
    803        
    804    
     802
     803
    805804    return a, b
    806        
    807                    
    808 
    809 def limit(quantity):       
     805
     806
     807
     808def limit(quantity):
    810809    """Limit slopes for each volume to eliminate artificial variance
    811810    introduced by e.g. second order extrapolator
    812    
     811
    813812    This is an unsophisticated limiter as it does not take into
    814813    account dependencies among quantities.
    815    
     814
    816815    precondition:
    817816    vertex values are estimated from gradient
    … …  
    823822
    824823    N = quantity.domain.number_of_elements
    825    
     824
    826825    beta_w = quantity.domain.beta_w
    827        
     826
    828827    qc = quantity.centroid_values
    829828    qv = quantity.vertex_values
    830        
     829
    831830    #Find min and max of this and neighbour's centroid values
    832831    qmax = zeros(qc.shape, Float)
    833     qmin = zeros(qc.shape, Float)       
    834    
     832    qmin = zeros(qc.shape, Float)
     833
    835834    for k in range(N):
    836835        qmax[k] = qmin[k] = qc[k]
    … …  
    839838            if n >= 0:
    840839                qn = qc[n] #Neighbour's centroid value
    841                
     840
    842841                qmin[k] = min(qmin[k], qn)
    843842                qmax[k] = max(qmax[k], qn)
    844      
    845    
     843
     844
    846845    #Diffences between centroids and maxima/minima
    847     dqmax = qmax - qc 
     846    dqmax = qmax - qc
    848847    dqmin = qmin - qc
    849        
     848
    850849    #Deltas between vertex and centroid values
    851850    dq = zeros(qv.shape, Float)
    … …  
    853852        dq[:,i] = qv[:,i] - qc
    854853
    855     #Phi limiter   
     854    #Phi limiter
    856855    for k in range(N):
    857        
     856
    858857        #Find the gradient limiter (phi) across vertices
    859858        phi = 1.0
    … …  
    862861            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
    863862            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
    864            
    865             phi = min( min(r*beta_w, 1), phi )   
     863
     864            phi = min( min(r*beta_w, 1), phi )
    866865
    867866        #Then update using phi limiter
    868         for i in range(3):           
     867        for i in range(3):
    869868            qv[k,i] = qc[k] + phi*dq[k,i]
    870869
    … …  
    874873if compile.can_use_C_extension('quantity_ext.c'):
    875874    #Replace python version with c implementations
    876        
     875
    877876    from quantity_ext import limit, compute_gradients,\
    878     extrapolate_second_order, interpolate_from_vertices_to_edges, update
    879 
     877    extrapolate_second_order, interpolate_from_vertices_to_edges, update
Note: See TracChangeset for help on using the changeset viewer.