source: development/pyvolution-1d/quantity.py @ 2716

Last change on this file since 2716 was 2716, checked in by jakeman, 18 years ago

Adding new files

File size: 20.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
18class Quantity:
19
20    def __init__(self, domain, vertex_values=None):
21
22        from domain import Domain
23        from Numeric import array, zeros, Float
24
25        msg = 'First argument in Quantity.__init__ '
26        msg += 'must be of class Domain (or a subclass thereof)'
27        assert isinstance(domain, Domain), msg
28
29        if vertex_values is None:
30            N = domain.number_of_elements
31            self.vertex_values = zeros((N, 2), Float)
32        else:
33            self.vertex_values = array(vertex_values, Float)
34
35            N, V = self.vertex_values.shape
36            assert V == 2,\
37                   'Two vertex values per element must be specified'
38
39
40            msg = 'Number of vertex values (%d) must be consistent with'\
41                  %N
42            msg += 'number of elements in specified domain (%d).'\
43                   %domain.number_of_elements
44
45            assert N == domain.number_of_elements, msg
46
47        self.domain = domain
48
49        #Allocate space for other quantities
50        self.centroid_values = zeros(N, Float)
51        #self.edge_values = zeros((N, 2), Float)
52        #edge values are values of the ends of each interval
53       
54        #does oe dimension need edge values???
55
56        #Intialise centroid values
57        self.interpolate()
58
59    #Methods for operator overloading
60    def __len__(self):
61        return self.centroid_values.shape[0]
62   
63    def interpolate(self):
64        """Compute interpolated values at centroid
65        Pre-condition: vertex_values have been set
66        """
67
68        N = self.vertex_values.shape[0]
69        for i in range(N):
70            v0 = self.vertex_values[i, 0]
71            v1 = self.vertex_values[i, 1]
72
73            self.centroid_values[i] = (v0 + v1)/2
74
75    def set_values(self, X, location='vertices'):
76        """Set values for quantity
77
78        X: Compatible list, Numeric array (see below), constant or function
79        location: Where values are to be stored.
80                  Permissible options are: vertices, centroid
81                  Default is "vertices"
82
83        In case of location == 'centroid' the dimension values must
84        be a list of a Numerical array of length N, N being the number
85        of elements in the mesh. Otherwise it must be of dimension Nx3
86
87        The values will be stored in elements following their
88        internal ordering.
89
90        If values are described a function, it will be evaluated at specified points
91
92        If selected location is vertices, values for centroid and edges
93        will be assigned interpolated values.
94        In any other case, only values for the specified locations
95        will be assigned and the others will be left undefined.
96        """
97
98        if location not in ['vertices', 'centroids']:#, 'edges']:
99            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
100            raise msg
101
102        if X is None:
103            msg = 'Given values are None'
104            raise msg
105
106        import types
107
108        if callable(X):
109            #Use function specific method
110            self.set_function_values(X, location)
111        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
112            if location == 'centroids':
113                self.centroid_values[:] = X
114            else:
115                self.vertex_values[:] = X
116
117        else:
118            #Use array specific method
119            self.set_array_values(X, location)
120
121        if location == 'vertices':
122            #Intialise centroid and edge values
123            self.interpolate()
124
125
126
127
128    def set_function_values(self, f, location='vertices'):
129        """Set values for quantity using specified function
130
131        f: x -> z Function where x and z are arrays
132        location: Where values are to be stored.
133                  Permissible options are: vertices, centroid
134                  Default is "vertices"
135        """
136
137        if location == 'centroids':
138            P = self.domain.centroids
139            self.set_values(f(P), location)
140        else:
141            #Vertices
142            P = self.domain.get_vertices()
143            for i in range(2):
144                self.vertex_values[:,i] = f(P[:,i])
145
146
147    def set_array_values(self, values, location='vertices'):
148        """Set values for quantity
149
150        values: Numeric array
151        location: Where values are to be stored.
152                  Permissible options are: vertices, centroid, edges
153                  Default is "vertices"
154
155        In case of location == 'centroid' the dimension values must
156        be a list of a Numerical array of length N, N being the number
157        of elements in the mesh. Otherwise it must be of dimension Nx2
158
159        The values will be stored in elements following their
160        internal ordering.
161
162        If selected location is vertices, values for centroid
163        will be assigned interpolated values.
164        In any other case, only values for the specified locations
165        will be assigned and the others will be left undefined.
166        """
167
168        from Numeric import array, Float
169
170        values = array(values).astype(Float)
171
172        N = self.centroid_values.shape[0]
173
174        msg = 'Number of values must match number of elements'
175        assert values.shape[0] == N, msg
176
177        if location == 'centroids':
178            assert len(values.shape) == 1, 'Values array must be 1d'
179            self.centroid_values = values
180        #elif location == 'edges':
181        #    assert len(values.shape) == 2, 'Values array must be 2d'
182        #    msg = 'Array must be N x 2'
183        #    self.edge_values = values
184        else:
185            assert len(values.shape) == 2, 'Values array must be 2d'
186            msg = 'Array must be N x 2'
187            assert values.shape[1] == 2, msg
188
189            self.vertex_values = values
190
191
192    def get_values(self, location='vertices', indices = None):
193        """get values for quantity
194
195        return X, Compatible list, Numeric array (see below)
196        location: Where values are to be stored.
197                  Permissible options are: vertices, edges, centroid
198                  and unique vertices. Default is 'vertices'
199
200        In case of location == 'centroids' the dimension values must
201        be a list of a Numerical array of length N, N being the number
202        of elements. Otherwise it must be of dimension Nx3
203
204        The returned values with be a list the length of indices
205        (N if indices = None).  Each value will be a list of the three
206        vertex values for this quantity.
207
208        Indices is the set of element ids that the operation applies to.
209
210        """
211        from Numeric import take
212
213        #if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
214        if location not in ['vertices', 'centroids', 'unique vertices']:
215            msg = 'Invalid location: %s' %location
216            raise msg
217
218        import types, Numeric
219        assert type(indices) in [types.ListType, types.NoneType,
220                                 Numeric.ArrayType],\
221                                 'Indices must be a list or None'
222
223        if location == 'centroids':
224            if (indices ==  None):
225                indices = range(len(self))
226            return take(self.centroid_values,indices)
227        #elif location == 'edges':
228        #    if (indices ==  None):
229        #        indices = range(len(self))
230        #    return take(self.edge_values,indices)
231        elif location == 'unique vertices':
232            if (indices ==  None):
233                indices=range(self.domain.coordinates.shape[0])
234            vert_values = []
235            #Go through list of unique vertices
236            for unique_vert_id in indices:
237                triangles = self.domain.vertexlist[unique_vert_id]
238
239                #In case there are unused points
240                if triangles is None:
241                    msg = 'Unique vertex not associated with triangles'
242                    raise msg
243
244                # Go through all triangle, vertex pairs
245                # Average the values
246                sum = 0
247                for triangle_id, vertex_id in triangles:
248                    sum += self.vertex_values[triangle_id, vertex_id]
249                vert_values.append(sum/len(triangles))
250            return Numeric.array(vert_values)
251        else:
252            if (indices ==  None):
253                indices = range(len(self))
254            return take(self.vertex_values,indices)
255
256    #Method for outputting model results
257    #FIXME: Split up into geometric and numeric stuff.
258    #FIXME: Geometric (X,Y,V) should live in mesh.py
259    #FIXME: STill remember to move XY to mesh
260    def get_vertex_values(self,
261                          #xy=True,
262                          x=True,
263                          smooth = None,
264                          precision = None,
265                          reduction = None):
266        """Return vertex values like an OBJ format
267
268        The vertex values are returned as one sequence in the 1D float array A.
269        If requested the coordinates will be returned in 1D arrays X and Y.
270
271        The connectivity is represented as an integer array, V, of dimension
272        M x 3, where M is the number of volumes. Each row has three indices
273        into the X, Y, A arrays defining the triangle.
274
275        if smooth is True, vertex values corresponding to one common
276        coordinate set will be smoothed according to the given
277        reduction operator. In this case vertex coordinates will be
278        de-duplicated.
279
280        If no smoothings is required, vertex coordinates and values will
281        be aggregated as a concatenation of values at
282        vertices 0, vertices 1 and vertices 2
283
284
285        Calling convention
286        if xy is True:
287           X,Y,A,V = get_vertex_values
288        else:
289           A,V = get_vertex_values
290
291        """
292
293        from Numeric import concatenate, zeros, Float, Int, array, reshape
294
295
296        if smooth is None:
297            smooth = self.domain.smooth
298
299        if precision is None:
300            precision = Float
301
302        if reduction is None:
303            reduction = self.domain.reduction
304
305        #Create connectivity
306
307        if smooth == True:
308
309            V = self.domain.get_vertices()
310            N = len(self.domain.vertexlist)
311            #N = len(self.domain.vertices)
312            A = zeros(N, precision)
313
314            #Smoothing loop
315            for k in range(N):
316                L = self.domain.vertexlist[k]
317                #L = self.domain.vertices[k]
318
319                #Go through all triangle, vertex pairs
320                #contributing to vertex k and register vertex value
321
322                if L is None: continue #In case there are unused points
323
324                contributions = []
325                for volume_id, vertex_id in L:
326                    v = self.vertex_values[volume_id, vertex_id]
327                    contributions.append(v)
328
329                A[k] = reduction(contributions)
330
331
332            #if xy is True:
333            if x is True:
334                 #X = self.domain.coordinates[:,0].astype(precision)
335                 X = self.domain.coordinates[:].astype(precision)
336                 #Y = self.domain.coordinates[:,1].astype(precision)
337
338                 #return X, Y, A, V
339                 return X, A, V
340           
341            #else:
342            return A, V
343        else:
344            #Don't smooth
345            #obj machinery moved to general_mesh
346
347            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
348            # These vert_id's will relate to the verts created below
349            #m = len(self.domain)  #Number of volumes
350            #M = 3*m        #Total number of unique vertices
351            #V = reshape(array(range(M)).astype(Int), (m,3))
352
353            #V = self.domain.get_triangles(obj=True)
354            V = self.domain.get_vertices
355            #FIXME use get_vertices, when ready
356
357            A = self.vertex_values.flat
358
359            #Do vertex coordinates
360            #if xy is True:
361            if x is True:
362                C = self.domain.get_vertex_coordinates()
363
364                X = C[:,0:6:2].copy()
365                Y = C[:,1:6:2].copy()
366
367                return X.flat, Y.flat, A, V
368            else:
369                return A, V
370
371
372class Conserved_quantity(Quantity):
373    """Class conserved quantity adds to Quantity:
374
375    storage and method for updating, and
376    methods for extrapolation from centropid to vertices inluding
377    gradients and limiters
378    """
379
380    def __init__(self, domain, vertex_values=None):
381        Quantity.__init__(self, domain, vertex_values)
382
383        from Numeric import zeros, Float
384       
385        #Allocate space for boundary values
386        #L = len(domain.boundary)
387        self.boundary_values = zeros(2, Float) #assumes no parrellism
388
389        #Allocate space for updates of conserved quantities by
390        #flux calculations and forcing functions
391       
392        N = domain.number_of_elements
393        self.explicit_update = zeros(N, Float )
394        self.semi_implicit_update = zeros(N, Float )
395
396        self.gradients = zeros(N, Float)
397        self.qmax = zeros(self.centroid_values.shape, Float)
398        self.qmin = zeros(self.centroid_values.shape, Float)
399
400        self.beta = domain.beta
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 sum, equal, ones, Float
409
410        N = self.centroid_values.shape[0]
411
412        #Explicit updates
413        self.centroid_values += timestep*self.explicit_update
414
415        #Semi implicit updates
416        denominator = ones(N, Float)-timestep*self.semi_implicit_update
417
418        if sum(equal(denominator, 0.0)) > 0.0:
419            msg = 'Zero division in semi implicit update. Call Stephen :-)'
420            raise msg
421        else:
422            #Update conserved_quantities from semi implicit updates
423            self.centroid_values /= denominator
424
425
426    def compute_gradients(self):
427        """Compute gradients of piecewise linear function defined by centroids of
428        neighbouring volumes.
429        """
430
431
432        from Numeric import array, zeros, Float
433
434        N = self.centroid_values.shape[0]
435
436
437        G = self.gradients
438        Q = self.centroid_values
439        X = self.domain.centroids
440
441        for k in range(N):
442
443            # first and last elements have boundaries
444
445            if k == 0:
446
447                #Get data
448                k0 = k
449                k1 = k+1
450
451                q0 = Q[k0]
452                q1 = Q[k1]
453
454                x0 = X[k0] #V0 centroid
455                x1 = X[k1] #V1 centroid
456
457                #Gradient
458                G[k] = (q1 - q0)/(x1 - x0)
459
460            elif k == N-1:
461
462                #Get data
463                k0 = k
464                k1 = k-1
465
466                q0 = Q[k0]
467                q1 = Q[k1]
468
469                x0 = X[k0] #V0 centroid
470                x1 = X[k1] #V1 centroid
471
472                #Gradient
473                G[k] = (q1 - q0)/(x1 - x0)
474
475            else:
476                #Interior Volume (2 neighbours)
477
478                #Get data
479                k0 = k-1
480                k2 = k+1
481
482                q0 = Q[k0]
483                q1 = Q[k]
484                q2 = Q[k2]
485
486                x0 = X[k0] #V0 centroid
487                x1 = X[k] #V1 centroid (Self)
488                x2 = X[k2] #V2 centroid
489
490                #Gradient
491                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
492
493        return
494
495    def extrapolate_first_order(self):
496        """Extrapolate conserved quantities from centroid to
497        vertices for each volume using
498        first order scheme.
499        """
500
501        qc = self.centroid_values
502        qv = self.vertex_values
503
504        for i in range(2):
505            qv[:,i] = qc
506
507
508    def extrapolate_second_order(self):
509        """Extrapolate conserved quantities from centroid to
510        vertices for each volume using
511        second order scheme.
512        """
513
514        self.compute_gradients()
515
516        G = self.gradients
517        V = self.domain.vertices
518        Qc = self.centroid_values
519        Qv = self.vertex_values
520
521        #Check each triangle
522        for k in range(self.domain.number_of_elements):
523            #Centroid coordinates
524            x = self.domain.centroids[k]
525
526            #vertex coordinates
527            x0, x1 = V[k,:]
528
529            #Extrapolate
530            Qv[k,0] = Qc[k] + G[k]*(x0-x)
531            Qv[k,1] = Qc[k] + G[k]*(x1-x)
532
533
534
535
536    def limit(self):
537        """Limit slopes for each volume to eliminate artificial variance
538        introduced by e.g. second order extrapolator
539
540        This is an unsophisticated limiter as it does not take into
541        account dependencies among quantities.
542
543        precondition:
544        vertex values are estimated from gradient
545        postcondition:
546        vertex values are updated
547        """
548
549        from Numeric import zeros, Float
550
551        N = self.domain.number_of_elements
552        beta = self.beta
553
554        qc = self.centroid_values
555        qv = self.vertex_values
556
557        #Find min and max of this and neighbour's centroid values
558        qmax = self.qmax
559        qmin = self.qmin
560
561        for k in range(N):
562            qmax[k] = qmin[k] = qc[k]
563
564            for i in [-1,1]:
565                n = k+i
566                if (n >= 0) & (n <= N-1):
567                    qn = qc[n] #Neighbour's centroid value
568
569                    qmin[k] = min(qmin[k], qn)
570                    qmax[k] = max(qmax[k], qn)
571
572        #Phi limiter
573        for k in range(N):
574
575            #Diffences between centroids and maxima/minima
576            dqmax = qmax[k] - qc[k]
577            dqmin = qmin[k] - qc[k]
578
579            #Deltas between vertex and centroid values
580            dq = [0.0, 0.0]
581            for i in range(2):
582                dq[i] = qv[k,i] - qc[k]
583
584            #Find the gradient limiter (phi) across vertices
585            phi = 1.0
586            for i in range(2):
587                r = 1.0
588                if (dq[i] > 0): r = dqmax/dq[i]
589                if (dq[i] < 0): r = dqmin/dq[i]
590
591                phi = min( min(r*beta, 1), phi )
592
593            #Then update using phi limiter
594            for i in range(2):
595                qv[k,i] = qc[k] + phi*dq[i]
596
597
598
599if __name__ == "__main__":
600    from domain import Domain
601
602    points1 = [0.0, 1.0, 2.0, 3.0]
603    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
604
605    D1 = Domain(points1)
606
607    Q1 = Conserved_quantity(D1, vertex_values)
608
609    print Q1.vertex_values
610    print Q1.centroid_values
611
612    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
613
614    Q1.set_values(new_vertex_values)
615
616    print Q1.vertex_values
617    print Q1.centroid_values
618
619    new_centroid_values = [20,30,40]
620    Q1.set_values(new_centroid_values,'centroids')
621
622    print Q1.vertex_values
623    print Q1.centroid_values
624
625    def fun(x):
626        return x**2
627
628    Q1.set_values(fun,'vertices')
629
630    print Q1.vertex_values
631    print Q1.centroid_values
632
633    Xc = Q1.domain.vertices
634    Qc = Q1.vertex_values
635    print Xc
636    print Qc
637
638    Qc[1,0] = 3
639
640    Q1.beta = 1.0
641    Q1.extrapolate_second_order()
642    Q1.limit()
643
644    from pylab import *
645    plot(Xc,Qc)
646
647    show()
648
649    points2 = arange(10)
650
651
652    D2 = Domain(points2)
653
654    Q2 = Conserved_quantity(D2)
655    Q2.set_values(fun,'vertices')
656    Xc = Q2.domain.vertices
657    Qc = Q2.vertex_values
658    plot(Xc,Qc)
659
660    show()
661    Q2.extrapolate_second_order()
662    Q2.limit()
663    Xc = Q2.domain.vertices
664    Qc = Q2.vertex_values
665
666    print Q2.centroid_values
667    print Qc
668
669    plot(Xc,Qc)
670
671    show()
672
673
674
675
676
677
678
679
Note: See TracBrowser for help on using the repository browser.