Changeset 5858


Ignore:
Timestamp:
Oct 22, 2008, 7:36:23 PM (15 years ago)
Author:
steve
Message:
 
Location:
anuga_work/development/anuga_1d
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/quantity.py

    r5844 r5858  
    1616
    1717
     18
    1819class Quantity:
    1920
     21   
    2022    def __init__(self, domain, vertex_values=None):
    21 
    22         #from domain import Domain
    23         #from domain_order2 import Domain
     23        #Initialise Quantity using optional vertex values.
     24       
    2425        from domain import Domain
    2526        from Numeric import array, zeros, Float
     
    7879        self.beta = domain.beta       
    7980
    80     #Methods for operator overloading
     81
    8182    def __len__(self):
     83        """
     84        Returns number of intervals.
     85        """
    8286        return self.centroid_values.shape[0]
    8387   
    8488    def interpolate(self):
    85         """Compute interpolated values at centroid
     89        """
     90        Compute interpolated values at centroid
    8691        Pre-condition: vertex_values have been set
    8792        """
     
    104109        In case of location == 'centroid' the dimension values must
    105110        be a list of a Numerical array of length N, N being the number
    106         of elements in the mesh. Otherwise it must be of dimension Nx3
     111        of elements in the mesh. Otherwise it must be of dimension Nx2
    107112
    108113        The values will be stored in elements following their
     
    117122        """
    118123
    119         if location not in ['vertices', 'centroids']:#, 'edges']:
     124        if location not in ['vertices', 'centroids']:
    120125            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
    121126            raise msg
     
    143148            #Intialise centroid and edge values
    144149            self.interpolate()
     150
    145151
    146152
     
    217223        return X, Compatible list, Numeric array (see below)
    218224        location: Where values are to be stored.
    219                   Permissible options are: vertices, edges, centroid
     225                  Permissible options are: vertices, centroid
    220226                  and unique vertices. Default is 'vertices'
    221227
     
    233239        from Numeric import take
    234240
    235         #if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
    236241        if location not in ['vertices', 'centroids', 'unique vertices']:
    237242            msg = 'Invalid location: %s' %location
     
    247252                indices = range(len(self))
    248253            return take(self.centroid_values,indices)
    249         #elif location == 'edges':
    250         #    if (indices ==  None):
    251         #        indices = range(len(self))
    252         #    return take(self.edge_values,indices)
    253254        elif location == 'unique vertices':
    254255            if (indices ==  None):
     
    257258            #Go through list of unique vertices
    258259            for unique_vert_id in indices:
    259                 triangles = self.domain.vertexlist[unique_vert_id]
     260                cells = self.domain.vertexlist[unique_vert_id]
    260261
    261262                #In case there are unused points
    262                 if triangles is None:
    263                     msg = 'Unique vertex not associated with triangles'
     263                if cells is None:
     264                    msg = 'Unique vertex not associated with cells'
    264265                    raise msg
    265266
    266                 # Go through all triangle, vertex pairs
     267                # Go through all cells, vertex pairs
    267268                # Average the values
    268269                sum = 0
    269                 for triangle_id, vertex_id in triangles:
    270                     sum += self.vertex_values[triangle_id, vertex_id]
    271                 vert_values.append(sum/len(triangles))
     270                for cell_id, vertex_id in cells:
     271                    sum += self.vertex_values[cell_id, vertex_id]
     272                vert_values.append(sum/len(cells))
    272273            return Numeric.array(vert_values)
    273274        else:
     
    276277            return take(self.vertex_values,indices)
    277278
    278     #Method for outputting model results
    279     #FIXME: Split up into geometric and numeric stuff.
    280     #FIXME: Geometric (X,Y,V) should live in mesh.py
    281     #FIXME: STill remember to move XY to mesh
     279
    282280    def get_vertex_values(self,
    283                           #xy=True,
    284281                          x=True,
    285282                          smooth = None,
     
    289286
    290287        The vertex values are returned as one sequence in the 1D float array A.
    291         If requested the coordinates will be returned in 1D arrays X and Y.
     288        If requested the coordinates will be returned in 1D arrays X.
    292289
    293290        The connectivity is represented as an integer array, V, of dimension
    294         M x 3, where M is the number of volumes. Each row has three indices
    295         into the X, Y, A arrays defining the triangle.
     291        M x 2, where M is the number of volumes. Each row has two indices
     292        into the X, A arrays defining the element.
    296293
    297294        if smooth is True, vertex values corresponding to one common
     
    302299        If no smoothings is required, vertex coordinates and values will
    303300        be aggregated as a concatenation of values at
    304         vertices 0, vertices 1 and vertices 2
     301        vertices 0, vertices 1
    305302
    306303
    307304        Calling convention
    308         if xy is True:
    309            X,Y,A,V = get_vertex_values
     305        if x is True:
     306           X,A,V = get_vertex_values
    310307        else:
    311308           A,V = get_vertex_values
     
    351348                A[k] = reduction(contributions)
    352349
    353 
    354             #if xy is True:
    355350            if x is True:
    356351                 #X = self.domain.coordinates[:,0].astype(precision)
     
    380375
    381376            #Do vertex coordinates
    382             #if xy is True:
    383377            if x is True:
    384                 C = self.domain.get_vertex_coordinates()
    385 
    386                 X = C[:,0:6:2].copy()
    387                 Y = C[:,1:6:2].copy()
    388 
    389                 return X.flat, Y.flat, A, V
     378                X = self.domain.get_vertex_coordinates()
     379
     380                #X = C[:,0:6:2].copy()
     381                #Y = C[:,1:6:2].copy()
     382
     383                return X.flat, A, V
    390384            else:
    391385                return A, V
     
    412406        N = self.centroid_values.shape[0]
    413407
    414         #print "pre update",self.centroid_values
    415         #print "timestep",timestep
    416         #print min(self.domain.areas)
    417         #print "explicit_update", self.explicit_update
    418 
    419408        #Explicit updates
    420409        self.centroid_values += timestep*self.explicit_update
    421 
    422         #print "post update",self.centroid_values
    423410       
    424411        #Semi implicit updates
     
    438425        """
    439426
     427        #print 'compute_gradient'
    440428
    441429        from Numeric import array, zeros, Float
     
    457445                k0 = k
    458446                k1 = k+1
     447                k2 = k+2
    459448
    460449                q0 = Q[k0]
    461450                q1 = Q[k1]
     451                q2 = Q[k2]
    462452
    463453                x0 = X[k0] #V0 centroid
    464454                x1 = X[k1] #V1 centroid
     455                x2 = X[k2]
    465456
    466457                #Gradient
    467                 G[k] = (q1 - q0)/(x1 - x0)
     458                #G[k] = (q1 - q0)/(x1 - x0)
     459               
     460                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
     461                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
    468462
    469463            elif k == N-1:
     
    472466                k0 = k
    473467                k1 = k-1
     468                k2 = k-2
    474469
    475470                q0 = Q[k0]
    476471                q1 = Q[k1]
     472                q2 = Q[k2]
    477473
    478474                x0 = X[k0] #V0 centroid
    479475                x1 = X[k1] #V1 centroid
     476                x2 = X[k2]
    480477
    481478                #Gradient
    482                 G[k] = (q1 - q0)/(x1 - x0)
     479                #G[k] = (q1 - q0)/(x1 - x0)
     480               
     481                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
     482                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
     483
     484##                q0 = Q[k0]
     485##                q1 = Q[k1]
     486##
     487##                x0 = X[k0] #V0 centroid
     488##                x1 = X[k1] #V1 centroid
     489##
     490##                #Gradient
     491##                G[k] = (q1 - q0)/(x1 - x0)
    483492
    484493            else:
     
    507516        """
    508517
     518        #print 'compute_minmod_gradients'
     519       
    509520        from Numeric import array, zeros, Float,sign
    510521       
     
    533544                k0 = k
    534545                k1 = k+1
     546                k2 = k+2
    535547
    536548                q0 = Q[k0]
    537549                q1 = Q[k1]
     550                q2 = Q[k2]
    538551
    539552                x0 = X[k0] #V0 centroid
    540553                x1 = X[k1] #V1 centroid
     554                x2 = X[k2]
    541555
    542556                #Gradient
    543                 G[k] = (q1 - q0)/(x1 - x0)
     557                #G[k] = (q1 - q0)/(x1 - x0)
     558               
     559                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
     560                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
    544561
    545562            elif k == N-1:
     
    548565                k0 = k
    549566                k1 = k-1
     567                k2 = k-2
    550568
    551569                q0 = Q[k0]
    552570                q1 = Q[k1]
     571                q2 = Q[k2]
    553572
    554573                x0 = X[k0] #V0 centroid
    555574                x1 = X[k1] #V1 centroid
     575                x2 = X[k2]
    556576
    557577                #Gradient
    558                 G[k] = (q1 - q0)/(x1 - x0)
     578                #G[k] = (q1 - q0)/(x1 - x0)
     579               
     580                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
     581                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
     582
     583##                #Get data
     584##                k0 = k
     585##                k1 = k-1
     586##
     587##                q0 = Q[k0]
     588##                q1 = Q[k1]
     589##
     590##                x0 = X[k0] #V0 centroid
     591##                x1 = X[k1] #V1 centroid
     592##
     593##                #Gradient
     594##                G[k] = (q1 - q0)/(x1 - x0)
    559595
    560596            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
     
    732768        #print limiter
    733769       
     770        #print 'limit_range'
    734771        N = self.domain.number_of_elements
    735772        qc = self.centroid_values
     
    866903        # Call correct module function
    867904        # (either from this module or C-extension)
    868         #saxpy_centroid_values(self,a,b)
    869905        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
    870906       
     
    883919
    884920
     921##
     922##def newLinePlot(title='Simple Plot'):
     923##    import Gnuplot
     924##    g = Gnuplot.Gnuplot()
     925##    g.title(title)
     926##    g('set data style linespoints')
     927##    g.xlabel('x')
     928##    g.ylabel('y')
     929##    return g
     930##
     931##def linePlot(g,x,y):
     932##    import Gnuplot
     933##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
    885934
    886935def newLinePlot(title='Simple Plot'):
    887     import Gnuplot
    888     g = Gnuplot.Gnuplot()
     936    import pylab as g
     937    g.ion()
     938    g.hold(False)
    889939    g.title(title)
    890     g('set data style linespoints')
    891940    g.xlabel('x')
    892941    g.ylabel('y')
    893     return g
    894 
    895 def linePlot(g,x,y):
    896     import Gnuplot
    897     g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
    898 
    899 
    900 
    901 
    902 
     942   
     943
     944def linePlot(x,y):
     945    import pylab as g
     946    g.plot(x.flat,y.flat)
     947
     948
     949def closePlots():
     950    import pylab as g
     951    g.close('all')
     952   
    903953if __name__ == "__main__":
    904954    #from domain import Domain
     
    9501000    Qc[1,0] = 3
    9511001
    952     Q1.beta = 1.0
    9531002    Q1.extrapolate_second_order()
    954     Q1.limit_minmod()
    955 
    956     g1 = newLinePlot('plot 1')
    957     linePlot(g1,Xc,Qc)
     1003    #Q1.limit_minmod()
     1004
     1005    newLinePlot('plots')
     1006    linePlot(Xc,Qc)
     1007    raw_input('press return')
    9581008
    9591009    points2 = arange(10)
     
    9641014    Xc = Q2.domain.vertices
    9651015    Qc = Q2.vertex_values
    966 
    967     g2 = newLinePlot('plot 2')
    968     linePlot(g2,Xc,Qc)
    969 
     1016    linePlot(Xc,Qc)
     1017    raw_input('press return')
    9701018
    9711019   
    9721020    Q2.extrapolate_second_order()
    973     Q2.limit_minmod()
     1021    #Q2.limit_minmod()
    9741022    Xc = Q2.domain.vertices
    9751023    Qc = Q2.vertex_values
    976 
    9771024    print Q2.centroid_values
    9781025    print Qc
    979     raw_input('press_return')
    980    
    981     g3 = newLinePlot('plot 3')
    982     linePlot(g3,Xc,Qc)
     1026    linePlot(Xc,Qc)
    9831027    raw_input('press return')
    9841028
    9851029
    9861030    for i in range(10):
     1031        import pylab as g
     1032        g.hold(True)
    9871033        fun = FunClass(i/10.0)
    988         Q2.set_values(fun,'vertices')
     1034        Q2.set_values(fun,'centroids')
     1035        Q2.extrapolate_second_order()
     1036        #Q2.limit_minmod()
    9891037        Qc = Q2.vertex_values
    990         linePlot(g3,Xc,Qc)
    991 
    992     raw_input('press return')
    993 
     1038        linePlot(Xc,Qc)
     1039        raw_input('press return')
     1040
     1041    raw_input('press return to quit')
     1042closePlots()
  • anuga_work/development/anuga_1d/test_quantity.py

    r5536 r5858  
    66
    77from quantity import *
    8 from anuga.config import epsilon
     8from config import epsilon
    99from Numeric import allclose, array, ones, Float
    10 
    11 from anuga.fit_interpolate.fit import fit_to_mesh
    12 #from anuga.pyvolution.least_squares import fit_to_mesh         
     10       
    1311from domain import Domain
    14 from anuga.geospatial_data.geospatial_data import Geospatial_data
    15 from anuga.coordinate_transforms.geo_reference import Geo_reference
    16 from anuga.utilities.polygon import *
     12
    1713
    1814#Aux for fit_interpolate.fit example
Note: See TracChangeset for help on using the changeset viewer.