# source:trunk/anuga_work/development/sudi/sw_1d/shock_detector_shv/quantity.py@7906

Last change on this file since 7906 was 7906, checked in by steve, 13 years ago

File size: 31.1 KB
Line
1"""Class Quantity - Implements values at each 1d element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 2 array of values at each vertex for each element.
10                  Default None
11
12   If vertex_values are None Create array of zeros compatible with domain.
13   Otherwise check that it is compatible with dimenions of domain.
14   Otherwise raise an exception
15"""
16
17
18
19class Quantity:
20
21
22    def __init__(self, domain, vertex_values=None):
23        #Initialise Quantity using optional vertex values.
24
25        from domain import Domain
26        from Numeric import Float
27        from numpy import array, zeros
28
29        msg = 'First argument in Quantity.__init__ '
30        msg += 'must be of class Domain (or a subclass thereof)'
31        assert isinstance(domain, Domain), msg
32
33        if vertex_values is None:
34            N = domain.number_of_elements
35            self.vertex_values = zeros((N, 2), Float)
36        else:
37            self.vertex_values = array(vertex_values, Float)
38
39            N, V = self.vertex_values.shape
40            assert V == 2,\
41                   'Two vertex values per element must be specified'
42
43
44            msg = 'Number of vertex values (%d) must be consistent with'\
45                  %N
46            msg += 'number of elements in specified domain (%d).'\
47                   %domain.number_of_elements
48
49            assert N == domain.number_of_elements, msg
50
51        self.domain = domain
52
53        #Allocate space for other quantities
54        self.centroid_values = zeros(N, Float)
55        self.centroid_backup_values = zeros(N, Float)
56        #self.edge_values = zeros((N, 2), Float)
57        #edge values are values of the ends of each interval
58
59        #Intialise centroid values
60        self.interpolate()
61
62
63        from Numeric import Float
64        from numpy import zeros
65
66        #Allocate space for boundary values
67        #L = len(domain.boundary)
68        self.boundary_values = zeros(2, Float) #assumes no parrellism
69
70        #Allocate space for updates of conserved quantities by
71        #flux calculations and forcing functions
72
73        N = domain.number_of_elements
74        self.explicit_update = zeros(N, Float )
75        self.semi_implicit_update = zeros(N, Float )
76
78        self.qmax = zeros(self.centroid_values.shape, Float)
79        self.qmin = zeros(self.centroid_values.shape, Float)
80
81        self.beta = domain.beta
82
83
84    def __len__(self):
85        """
86        Returns number of intervals.
87        """
88        return self.centroid_values.shape[0]
89
90    def interpolate(self):
91        """
92        Compute interpolated values at centroid
93        Pre-condition: vertex_values have been set
94        """
95
96        N = self.vertex_values.shape[0]
97        for i in range(N):
98            v0 = self.vertex_values[i, 0]
99            v1 = self.vertex_values[i, 1]
100
101            self.centroid_values[i] = (v0 + v1)/2.0
102
103    def set_values(self, X, location='vertices'):
104        """Set values for quantity
105
106        X: Compatible list, Numeric array (see below), constant or function
107        location: Where values are to be stored.
108                  Permissible options are: vertices, centroid
109                  Default is "vertices"
110
111        In case of location == 'centroid' the dimension values must
112        be a list of a Numerical array of length N, N being the number
113        of elements in the mesh. Otherwise it must be of dimension Nx2
114
115        The values will be stored in elements following their
116        internal ordering.
117
118        If values are described a function, it will be evaluated at specified points
119
120        If selected location is vertices, values for centroid and edges
121        will be assigned interpolated values.
122        In any other case, only values for the specified locations
123        will be assigned and the others will be left undefined.
124        """
125
126        if location not in ['vertices', 'centroids']:
127            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
128            raise msg
129
130        if X is None:
131            msg = 'Given values are None'
132            raise msg
133
134        import types
135
136        if callable(X):
137            #Use function specific method
138            self.set_function_values(X, location)
139        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
140            if location == 'centroids':
141                self.centroid_values[:] = X
142            else:
143                self.vertex_values[:] = X
144
145        else:
146            #Use array specific method
147            self.set_array_values(X, location)
148
149        if location == 'vertices':
150            #Intialise centroid and edge values
151            self.interpolate()
152
153
154
155
156
157    def set_function_values(self, f, location='vertices'):
158        """Set values for quantity using specified function
159
160        f: x -> z Function where x and z are arrays
161        location: Where values are to be stored.
162                  Permissible options are: vertices, centroid
163                  Default is "vertices"
164        """
165
166        if location == 'centroids':
167
168            P = self.domain.centroids
169            self.set_values(f(P), location)
170        else:
171            #Vertices
172            P = self.domain.get_vertices()
173
174            for i in range(2):
175                self.vertex_values[:,i] = f(P[:,i])
176
177    def set_array_values(self, values, location='vertices'):
178        """Set values for quantity
179
180        values: Numeric array
181        location: Where values are to be stored.
182                  Permissible options are: vertices, centroid, edges
183                  Default is "vertices"
184
185        In case of location == 'centroid' the dimension values must
186        be a list of a Numerical array of length N, N being the number
187        of elements in the mesh. Otherwise it must be of dimension Nx2
188
189        The values will be stored in elements following their
190        internal ordering.
191
192        If selected location is vertices, values for centroid
193        will be assigned interpolated values.
194        In any other case, only values for the specified locations
195        will be assigned and the others will be left undefined.
196        """
197
198        from Numeric import Float
199        from numpy import array
200
201        values = array(values).astype(Float)
202
203        N = self.centroid_values.shape[0]
204
205        msg = 'Number of values must match number of elements'
206        assert values.shape[0] == N, msg
207
208        if location == 'centroids':
209            assert len(values.shape) == 1, 'Values array must be 1d'
210            self.centroid_values = values
211        #elif location == 'edges':
212        #    assert len(values.shape) == 2, 'Values array must be 2d'
213        #    msg = 'Array must be N x 2'
214        #    self.edge_values = values
215        else:
216            assert len(values.shape) == 2, 'Values array must be 2d'
217            msg = 'Array must be N x 2'
218            assert values.shape[1] == 2, msg
219
220            self.vertex_values = values
221
222
223    def get_values(self, location='vertices', indices = None):
224        """get values for quantity
225
226        return X, Compatible list, Numeric array (see below)
227        location: Where values are to be stored.
228                  Permissible options are: vertices, centroid
229                  and unique vertices. Default is 'vertices'
230
231        In case of location == 'centroids' the dimension values must
232        be a list of a Numerical array of length N, N being the number
233        of elements. Otherwise it must be of dimension Nx3
234
235        The returned values with be a list the length of indices
236        (N if indices = None).  Each value will be a list of the three
237        vertex values for this quantity.
238
239        Indices is the set of element ids that the operation applies to.
240
241        """
242        from Numeric import take
243
244        if location not in ['vertices', 'centroids', 'unique vertices']:
245            msg = 'Invalid location: %s' %location
246            raise msg
247
248        import types, Numeric
249        assert type(indices) in [types.ListType, types.NoneType,
250                                 Numeric.ArrayType],\
251                                 'Indices must be a list or None'
252
253        if location == 'centroids':
254            if (indices ==  None):
255                indices = range(len(self))
256            return take(self.centroid_values,indices)
257        elif location == 'unique vertices':
258            if (indices ==  None):
259                indices=range(self.domain.coordinates.shape[0])
260            vert_values = []
261            #Go through list of unique vertices
262            for unique_vert_id in indices:
263                cells = self.domain.vertexlist[unique_vert_id]
264
265                #In case there are unused points
266                if cells is None:
267                    msg = 'Unique vertex not associated with cells'
268                    raise msg
269
270                # Go through all cells, vertex pairs
271                # Average the values
272                sum = 0
273                for cell_id, vertex_id in cells:
274                    sum += self.vertex_values[cell_id, vertex_id]
275                vert_values.append(sum/len(cells))
276            return Numeric.array(vert_values)
277        else:
278            if (indices ==  None):
279                indices = range(len(self))
280            return take(self.vertex_values,indices)
281
282
283    def get_vertex_values(self,
284                          x=True,
285                          smooth = None,
286                          precision = None,
287                          reduction = None):
288        """Return vertex values like an OBJ format
289
290        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.
292
293        The connectivity is represented as an integer array, V, of dimension
294        M x 2, where M is the number of volumes. Each row has two indices
295        into the X, A arrays defining the element.
296
297        if smooth is True, vertex values corresponding to one common
298        coordinate set will be smoothed according to the given
299        reduction operator. In this case vertex coordinates will be
300        de-duplicated.
301
302        If no smoothings is required, vertex coordinates and values will
303        be aggregated as a concatenation of values at
304        vertices 0, vertices 1
305
306
307        Calling convention
308        if x is True:
309           X,A,V = get_vertex_values
310        else:
311           A,V = get_vertex_values
312
313        """
314
315        from Numeric import Float, Int
316        from numpy import concatenate, zeros, array, reshape
317
318
319        if smooth is None:
320            smooth = self.domain.smooth
321
322        if precision is None:
323            precision = Float
324
325        if reduction is None:
326            reduction = self.domain.reduction
327
328        #Create connectivity
329
330        if smooth == True:
331
332            V = self.domain.get_vertices()
333            N = len(self.domain.vertexlist)
334            #N = len(self.domain.vertices)
335            A = zeros(N, precision)
336
337            #Smoothing loop
338            for k in range(N):
339                L = self.domain.vertexlist[k]
340                #L = self.domain.vertices[k]
341
342                #Go through all triangle, vertex pairs
343                #contributing to vertex k and register vertex value
344
345                if L is None: continue #In case there are unused points
346
347                contributions = []
348                for volume_id, vertex_id in L:
349                    v = self.vertex_values[volume_id, vertex_id]
350                    contributions.append(v)
351
352                A[k] = reduction(contributions)
353
354            if x is True:
355                 #X = self.domain.coordinates[:,0].astype(precision)
356                 X = self.domain.coordinates[:].astype(precision)
357                 #Y = self.domain.coordinates[:,1].astype(precision)
358
359                 #return X, Y, A, V
360                 return X, A, V
361
362            #else:
363            return A, V
364        else:
365            #Don't smooth
366            #obj machinery moved to general_mesh
367
368            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
369            # These vert_id's will relate to the verts created below
370            #m = len(self.domain)  #Number of volumes
371            #M = 3*m        #Total number of unique vertices
372            #V = reshape(array(range(M)).astype(Int), (m,3))
373
374            #V = self.domain.get_triangles(obj=True)
375            V = self.domain.get_vertices
376            #FIXME use get_vertices, when ready
377
378            A = self.vertex_values.flat
379
380            #Do vertex coordinates
381            if x is True:
382                X = self.domain.get_vertex_coordinates()
383
384                #X = C[:,0:6:2].copy()
385                #Y = C[:,1:6:2].copy()
386
387                return X.flat, A, V
388            else:
389                return A, V
390
391    def get_integral(self):
392        """Compute the integral of quantity across entire domain
393        """
394        integral = 0
395        for k in range(self.domain.number_of_elements):
396            area = self.domain.areas[k]
397            qc = self.centroid_values[k]
398            integral += qc*area
399
400        return integral
401
402
403    def update(self, timestep):
404        """Update centroid values based on values stored in
405        explicit_update and semi_implicit_update as well as given timestep
406        """
407
408        from Numeric import Float
409        from numpy import sum, equal, ones
410
411        N = self.centroid_values.shape[0]
412
414        self.centroid_values += timestep*self.explicit_update
415
417        denominator = ones(N, Float)-timestep*self.semi_implicit_update
418
419        if sum(equal(denominator, 0.0)) > 0.0:
420            msg = 'Zero division in semi implicit update. Call Stephen :-)'
421            raise msg
422        else:
423            #Update conserved_quantities from semi implicit updates
424            self.centroid_values /= denominator
425
426
428        """Compute gradients of piecewise linear function defined by centroids of
429        neighbouring volumes.
430        """
431
433
434        from Numeric import Float
435        from numpy import array, zeros
436
437        N = self.centroid_values.shape[0]
438
439
441        Q = self.centroid_values
442        X = self.domain.centroids
443
444        for k in range(N):
445
446            # first and last elements have boundaries
447
448            if k == 0:
449
450                #Get data
451                k0 = k
452                k1 = k+1
453                k2 = k+2
454
455                q0 = Q[k0]
456                q1 = Q[k1]
457                q2 = Q[k2]
458
459                x0 = X[k0] #V0 centroid
460                x1 = X[k1] #V1 centroid
461                x2 = X[k2]
462
464                #G[k] = (q1 - q0)/(x1 - x0)
465
466                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
467                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
468
469            elif k == N-1:
470
471                #Get data
472                k0 = k
473                k1 = k-1
474                k2 = k-2
475
476                q0 = Q[k0]
477                q1 = Q[k1]
478                q2 = Q[k2]
479
480                x0 = X[k0] #V0 centroid
481                x1 = X[k1] #V1 centroid
482                x2 = X[k2]
483
485                #G[k] = (q1 - q0)/(x1 - x0)
486
487                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
488                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
489
490##                q0 = Q[k0]
491##                q1 = Q[k1]
492##
493##                x0 = X[k0] #V0 centroid
494##                x1 = X[k1] #V1 centroid
495##
497##                G[k] = (q1 - q0)/(x1 - x0)
498
499            else:
500                #Interior Volume (2 neighbours)
501
502                #Get data
503                k0 = k-1
504                k2 = k+1
505
506                q0 = Q[k0]
507                q1 = Q[k]
508                q2 = Q[k2]
509
510                x0 = X[k0] #V0 centroid
511                x1 = X[k]  #V1 centroid (Self)
512                x2 = X[k2] #V2 centroid
513
515                #G[k] = (q2-q0)/(x2-x0)
516                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
517
518
520        """Compute gradients of piecewise linear function defined by centroids of
521        neighbouring volumes.
522        """
523
525
526        from Numeric import Float,sign
527        from numpy import array, zeros
528
529        def xmin(a,b):
530            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
531
532        def xmic(t,a,b):
533            return xmin(t*xmin(a,b), 0.50*(a+b) )
534
535
536
537        N = self.centroid_values.shape[0]
538
539
541        Q = self.centroid_values
542        X = self.domain.centroids
543
544        for k in range(N):
545
546            # first and last elements have boundaries
547
548            if k == 0:
549
550                #Get data
551                k0 = k
552                k1 = k+1
553                k2 = k+2
554
555                q0 = Q[k0]
556                q1 = Q[k1]
557                q2 = Q[k2]
558
559                x0 = X[k0] #V0 centroid
560                x1 = X[k1] #V1 centroid
561                x2 = X[k2]
562
564                #G[k] = (q1 - q0)/(x1 - x0)
565
566                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
567                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
568
569            elif k == N-1:
570
571                #Get data
572                k0 = k
573                k1 = k-1
574                k2 = k-2
575
576                q0 = Q[k0]
577                q1 = Q[k1]
578                q2 = Q[k2]
579
580                x0 = X[k0] #V0 centroid
581                x1 = X[k1] #V1 centroid
582                x2 = X[k2]
583
585                #G[k] = (q1 - q0)/(x1 - x0)
586
587                G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0)
588                G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1)
589
590##                #Get data
591##                k0 = k
592##                k1 = k-1
593##
594##                q0 = Q[k0]
595##                q1 = Q[k1]
596##
597##                x0 = X[k0] #V0 centroid
598##                x1 = X[k1] #V1 centroid
599##
601##                G[k] = (q1 - q0)/(x1 - x0)
602
603            elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
604                G[k] = 0.0
605
606            else:
607                #Interior Volume (2 neighbours)
608
609                #Get data
610                k0 = k-1
611                k2 = k+1
612
613                q0 = Q[k0]
614                q1 = Q[k]
615                q2 = Q[k2]
616
617                x0 = X[k0] #V0 centroid
618                x1 = X[k]  #V1 centroid (Self)
619                x2 = X[k2] #V2 centroid
620
621                # assuming uniform grid
622                d1 = (q1 - q0)/(x1-x0)
623                d2 = (q2 - q1)/(x2-x1)
624
626                #G[k] = (d1+d2)*0.5
627                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)
628                G[k] = xmic( self.domain.beta, d1, d2 )
629
630
631    def extrapolate_first_order(self):
632        """Extrapolate conserved quantities from centroid to
633        vertices for each volume using
634        first order scheme.
635        """
636
637        qc = self.centroid_values
638        qv = self.vertex_values
639
640        for i in range(2):
641            qv[:,i] = qc
642
643
644    def extrapolate_second_order(self):
645        """Extrapolate conserved quantities from centroid to
646        vertices for each volume using
647        second order scheme.
648        """
649        if self.domain.limiter == "pyvolution":
654
659
661            V = self.domain.vertices
662            qc = self.centroid_values
663            qv = self.vertex_values
664
665            #Check each triangle
666            for k in range(self.domain.number_of_elements):
667                #Centroid coordinates
668                x = self.domain.centroids[k]
669
670                #vertex coordinates
671                x0, x1 = V[k,:]
672
673                #Extrapolate
674                qv[k,0] = qc[k] + G[k]*(x0-x)
675                qv[k,1] = qc[k] + G[k]*(x1-x)
676            self.limit_pyvolution()
677        elif self.domain.limiter == "minmod_steve":
678            self.limit_minmod()
679        else:
680            self.limit_range()
681
682
683
684    def limit_minmod(self):
689
691        V = self.domain.vertices
692        qc = self.centroid_values
693        qv = self.vertex_values
694
695        #Check each triangle
696        for k in range(self.domain.number_of_elements):
697            #Centroid coordinates
698            x = self.domain.centroids[k]
699
700            #vertex coordinates
701            x0, x1 = V[k,:]
702
703            #Extrapolate
704            qv[k,0] = qc[k] + G[k]*(x0-x)
705            qv[k,1] = qc[k] + G[k]*(x1-x)
706
707
708    def limit_pyvolution(self):
709        """
710        Limit slopes for each volume to eliminate artificial variance
711        introduced by e.g. second order extrapolator
712
713        This is an unsophisticated limiter as it does not take into
714        account dependencies among quantities.
715
716        precondition:
717        vertex values are estimated from gradient
718        postcondition:
719        vertex values are updated
720        """
721        from Numeric import Float
722        from numpy import zeros
723
724        N = self.domain.number_of_elements
725        beta = self.domain.beta
726        #beta = 0.8
727
728        qc = self.centroid_values
729        qv = self.vertex_values
730
731        #Find min and max of this and neighbour's centroid values
732        qmax = self.qmax
733        qmin = self.qmin
734
735        for k in range(N):
736            qmax[k] = qmin[k] = qc[k]
737            for i in range(2):
738                n = self.domain.neighbours[k,i]
739                if n >= 0:
740                    qn = qc[n] #Neighbour's centroid value
741
742                    qmin[k] = min(qmin[k], qn)
743                    qmax[k] = max(qmax[k], qn)
744
745
746        #Diffences between centroids and maxima/minima
747        dqmax = qmax - qc
748        dqmin = qmin - qc
749
750        #Deltas between vertex and centroid values
751        dq = zeros(qv.shape, Float)
752        for i in range(2):
753            dq[:,i] = qv[:,i] - qc
754
755        #Phi limiter
756        for k in range(N):
757
758            #Find the gradient limiter (phi) across vertices
759            phi = 1.0
760            for i in range(2):
761                r = 1.0
762                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
763                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
764
765                phi = min( min(r*beta, 1), phi )
766
767            #Then update using phi limiter
768            for i in range(2):
769                qv[k,i] = qc[k] + phi*dq[k,i]
770
771    def limit_range(self):
772        import sys
773        from Numeric import Float
774        from numpy import zeros
775        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
776        limiter = self.domain.limiter
777        #print limiter
778
779        #print 'limit_range'
780        N = self.domain.number_of_elements
781        qc = self.centroid_values
782        qv = self.vertex_values
783        C = self.domain.centroids
784        X = self.domain.vertices
785        beta_p = zeros(N,Float)
786        beta_m = zeros(N,Float)
787        beta_x = zeros(N,Float)
788
789        for k in range(N):
790
791            n0 = self.domain.neighbours[k,0]
792            n1 = self.domain.neighbours[k,1]
793
794            if ( n0 >= 0) & (n1 >= 0):
795                #SLOPE DERIVATIVE LIMIT
796                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
797                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
798                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
799
800        dq = zeros(qv.shape, Float)
801        for i in range(2):
802            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
803
804        #Phi limiter
805        for k in range(N):
806            n0 = self.domain.neighbours[k,0]
807            n1 = self.domain.neighbours[k,1]
808            if n0 < 0:
809                #phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])    #This is the original
810                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
811            elif n1 < 0:
812                #phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])    #This is the original
813                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
814            else:
815                if limiter == "minmod":
816                    phi = minmod(beta_p[k],beta_m[k])  #(beta_p[k]+1.0e-15,beta_m[k])
817                elif limiter == "minmod_kurganov":#Change this
818                    # Also known as monotonized central difference limiter if theta = 2.0
819                    theta = 2.0
820                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
821                elif limiter == "superbee":
822                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
823                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
824                    phi = maxmod(slope1,slope2)
825                elif limiter == "vanleer":
826                    phi = vanleer(beta_p[k],beta_m[k])
829
830            for i in range(2):
831                qv[k,i] = qc[k] + phi*dq[k,i]
832
833    def limit_steve_slope(self):
834
835        import sys
836        from Numeric import Float
837        from numpy import zeros
838        from util import minmod, minmod_kurganov, maxmod, vanleer
839
840        N = self.domain.number_of_elements
841        limiter = self.domain.limiter
842        limiter_type = self.domain.limiter_type
843
844        qc = self.centroid_values
845        qv = self.vertex_values
846
847        #Find min and max of this and neighbour's centroid values
848        beta_p = zeros(N,Float)
849        beta_m = zeros(N,Float)
850        beta_x = zeros(N,Float)
851        C = self.domain.centroids
852        X = self.domain.vertices
853
854        for k in range(N):
855
856            n0 = self.domain.neighbours[k,0]
857            n1 = self.domain.neighbours[k,1]
858
859            if (n0 >= 0) & (n1 >= 0):
860                # Check denominator not zero
861                if (qc[k+1]-qc[k]) == 0.0:
862                    beta_p[k] = float(sys.maxint)
863                    beta_m[k] = float(sys.maxint)
864                else:
865                    #STEVE LIMIT
866                    beta_p[k] = (qc[k]-qc[k-1])/(qc[k+1]-qc[k])
867                    beta_m[k] = (qc[k+2]-qc[k+1])/(qc[k+1]-qc[k])
868
869        #Deltas between vertex and centroid values
870        dq = zeros(qv.shape, Float)
871        for i in range(2):
872            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
873
874        #Phi limiter
875        for k in range(N):
876
877            phi = 0.0
878            if limiter == "flux_minmod":
879                #FLUX MINMOD
880                phi = minmod_kurganov(1.0,beta_m[k],beta_p[k])
881            elif limiter == "flux_superbee":
882                #FLUX SUPERBEE
883                phi = max(0.0,min(1.0,2.0*beta_m[k]),min(2.0,beta_m[k]))+max(0.0,min(1.0,2.0*beta_p[k]),min(2.0,beta_p[k]))-1.0
884            elif limiter == "flux_muscl":
885                #FLUX MUSCL
886                phi = max(0.0,min(2.0,2.0*beta_m[k],2.0*beta_p[k],0.5*(beta_m[k]+beta_p[k])))
887            elif limiter == "flux_vanleer":
888                #FLUX VAN LEER
889                phi = (beta_m[k]+abs(beta_m[k]))/(1.0+abs(beta_m[k]))+(beta_p[k]+abs(beta_p[k]))/(1.0+abs(beta_p[k]))-1.0
890
891                #Then update using phi limiter
892                n = self.domain.neighbours[k,1]
893                if n>=0:
894                    #qv[k,0] = qc[k] - 0.5*phi*(qc[k+1]-qc[k])
895                    #qv[k,1] = qc[k] + 0.5*phi*(qc[k+1]-qc[k])
896                    qv[k,0] = qc[k] + 0.5*phi*(qv[k,0]-qc[k])
897                    qv[k,1] = qc[k] + 0.5*phi*(qv[k,1]-qc[k])
898                else:
899                    qv[k,i] = qc[k]
900
901    def backup_centroid_values(self):
902        # Call correct module function
903        # (either from this module or C-extension)
904        #backup_centroid_values(self)
905
906        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
907
908    def saxpy_centroid_values(self,a,b):
909        # Call correct module function
910        # (either from this module or C-extension)
911        self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
912
913class Conserved_quantity(Quantity):
914    """Class conserved quantity adds to Quantity:
915
916    storage and method for updating, and
917    methods for extrapolation from centropid to vertices inluding
919    """
920
921    def __init__(self, domain, vertex_values=None):
922        Quantity.__init__(self, domain, vertex_values)
923
924        print "Use Quantity instead of Conserved_quantity"
925
926"""
927## closePlots crashes with non-interactive matplotlib. (Sudi, 23 July 2010)
928##def newLinePlot(title='Simple Plot'):
929##    import Gnuplot
930##    g = Gnuplot.Gnuplot()
931##    g.title(title)
932##    g('set data style linespoints')
933##    g.xlabel('x')
934##    g.ylabel('y')
935##    return g
936##
937##def linePlot(g,x,y):
938##    import Gnuplot
939##    g.plot(Gnuplot.PlotItems.Data(x.flat,y.flat))
940
941def newLinePlot(title='Simple Plot'):
942    import pylab as g
943    g.ion()
944    g.hold(False)
945    g.title(title)
946    g.xlabel('x')
947    g.ylabel('y')
948
949
950def linePlot(x,y):
951    import pylab as g
952    g.plot(x.flat,y.flat)
953
954
955def closePlots():
956    import pylab as g
957    g.close('all')
958
959if __name__ == "__main__":
960    #from domain import Domain
961    from shallow_water_domain import Domain
962    from numpy import arange
963
964    points1 = [0.0, 1.0, 2.0, 3.0]
965    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
966
967    D1 = Domain(points1)
968
969    Q1 = Quantity(D1, vertex_values)
970
971    print Q1.vertex_values
972    print Q1.centroid_values
973
974    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
975
976    Q1.set_values(new_vertex_values)
977
978    print Q1.vertex_values
979    print Q1.centroid_values
980
981    new_centroid_values = [20,30,40]
982    Q1.set_values(new_centroid_values,'centroids')
983
984    print Q1.vertex_values
985    print Q1.centroid_values
986
987    class FunClass:
988        def __init__(self,value):
989            self.value = value
990
991        def __call__(self,x):
992            return self.value*(x**2)
993
994
995    fun = FunClass(1.0)
996    Q1.set_values(fun,'vertices')
997
998    print Q1.vertex_values
999    print Q1.centroid_values
1000
1001    Xc = Q1.domain.vertices
1002    Qc = Q1.vertex_values
1003    print Xc
1004    print Qc
1005
1006    Qc[1,0] = 3
1007
1008    Q1.extrapolate_second_order()
1009    #Q1.limit_minmod()
1010
1011    newLinePlot('plots')
1012    linePlot(Xc,Qc)
1013    raw_input('press return')
1014
1015    points2 = arange(10)
1016    D2 = Domain(points2)
1017
1018    Q2 = Quantity(D2)
1019    Q2.set_values(fun,'vertices')
1020    Xc = Q2.domain.vertices
1021    Qc = Q2.vertex_values
1022    linePlot(Xc,Qc)
1023    raw_input('press return')
1024
1025
1026    Q2.extrapolate_second_order()
1027    #Q2.limit_minmod()
1028    Xc = Q2.domain.vertices
1029    Qc = Q2.vertex_values
1030    print Q2.centroid_values
1031    print Qc
1032    linePlot(Xc,Qc)
1033    raw_input('press return')
1034
1035
1036    for i in range(10):
1037        import pylab as g
1038        g.hold(True)
1039        fun = FunClass(i/10.0)
1040        Q2.set_values(fun,'centroids')
1041        Q2.extrapolate_second_order()
1042        #Q2.limit_minmod()
1043        Qc = Q2.vertex_values
1044        linePlot(Xc,Qc)
1045        raw_input('press return')
1046