source: inundation/ga/storm_surge/pyvolution/quantity.py @ 1632

Last change on this file since 1632 was 1632, checked in by ole, 20 years ago

Interpolation work and some refactoring

File size: 29.2 KB
RevLine 
[1290]1"""Class Quantity - Implements values at each triangular element
2
3To create:
4
5   Quantity(domain, vertex_values)
6
7   domain: Associated domain structure. Required.
8
9   vertex_values: N x 3 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 mesh import Mesh
23        from Numeric import array, zeros, Float
24
25        msg = 'First argument in Quantity.__init__ '
26        msg += 'must be of class Mesh (or a subclass thereof)'
27        assert isinstance(domain, Mesh), msg
28
29        if vertex_values is None:
30            N = domain.number_of_elements
31            self.vertex_values = zeros((N, 3), Float)
32        else:
33            self.vertex_values = array(vertex_values).astype(Float)
34
35            N, V = self.vertex_values.shape
36            assert V == 3,\
37                   'Three 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, 3), Float)
52
53        #Intialise centroid and edge_values
54        self.interpolate()
55
56    def __len__(self):
57        return self.centroid_values.shape[0]
58
59    def interpolate(self):
60        """Compute interpolated values at edges and centroid
61        Pre-condition: vertex_values have been set
62        """
63
64        N = self.vertex_values.shape[0]
65        for i in range(N):
66            v0 = self.vertex_values[i, 0]
67            v1 = self.vertex_values[i, 1]
68            v2 = self.vertex_values[i, 2]
69
70            self.centroid_values[i] = (v0 + v1 + v2)/3
71
72        self.interpolate_from_vertices_to_edges()
73
74
75    def interpolate_from_vertices_to_edges(self):
76        #Call correct module function
77        #(either from this module or C-extension)
78        interpolate_from_vertices_to_edges(self)
79
80
81    def set_values(self, X, location='vertices', indexes = None):
82        """Set values for quantity
83
84        X: Compatible list, Numeric array (see below), constant or function
85        location: Where values are to be stored.
86                  Permissible options are: vertices, edges, centroids
87                  Default is "vertices"
88
89        In case of location == 'centroids' the dimension values must
90        be a list of a Numerical array of length N, N being the number
91        of elements. Otherwise it must be of dimension Nx3
92
93        The values will be stored in elements following their
94        internal ordering.
95
96        If values are described a function, it will be evaluated at
97        specified points
98
99        If indexex is not 'unique vertices' Indexes is the set of element ids
100        that the operation applies to.
101        If indexex is 'unique vertices' Indexes is the set of vertex ids
102        that the operation applies to.
103
104
105        If selected location is vertices, values for centroid and edges
106        will be assigned interpolated values.
107        In any other case, only values for the specified locations
108        will be assigned and the others will be left undefined.
109        """
110
111        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
112            msg = 'Invalid location: %s' %location
113            raise msg
114
115        if X is None:
116            msg = 'Given values are None'
117            raise msg
118
119        import types, Numeric
120        assert type(indexes) in [types.ListType, types.NoneType,
121                                 Numeric.ArrayType],\
122                                 'Indices must be a list or None'
123
124
125        if callable(X):
126            #Use function specific method
127            self.set_function_values(X, location, indexes = indexes)
128        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
129            if location == 'centroids':
130                if (indexes ==  None):
131                    self.centroid_values[:] = X
132                else:
133                    #Brute force
134                    for i in indexes:
135                        self.centroid_values[i,:] = X
136
137            elif location == 'edges':
138                if (indexes ==  None):
139                    self.edge_values[:] = X
140                else:
141                    #Brute force
142                    for i in indexes:
143                        self.edge_values[i,:] = X
144
145            elif location == 'unique vertices':
146                if (indexes ==  None):
147                    self.edge_values[:] = X
148                else:
149
150                    #Go through list of unique vertices
151                    for unique_vert_id in indexes:
152                        triangles = self.domain.vertexlist[unique_vert_id]
153
154                        #In case there are unused points
155                        if triangles is None: continue
156
157                        #Go through all triangle, vertex pairs
158                        #and set corresponding vertex value
159                        for triangle_id, vertex_id in triangles:
160                            self.vertex_values[triangle_id, vertex_id] = X
161
162                        #Intialise centroid and edge_values
163                        self.interpolate()
164            else:
165                if (indexes ==  None):
166                    self.vertex_values[:] = X
167                else:
168                    #Brute force
169                    for i_vertex in indexes:
170                        self.vertex_values[i_vertex,:] = X
171
172        else:
173            #Use array specific method
174            self.set_array_values(X, location, indexes = indexes)
175
176        if location == 'vertices' or location == 'unique vertices':
177            #Intialise centroid and edge_values
178            self.interpolate()
179
180        if location == 'centroids':
181            #Extrapolate 1st order - to capture notion of area being specified
182            self.extrapolate_first_order()
183
184    def get_values(self, location='vertices', indexes = None):
185        """get values for quantity
186
187        return X, Compatible list, Numeric array (see below)
188        location: Where values are to be stored.
189                  Permissible options are: vertices, edges, centroid
190                  Default is "vertices"
191
192        In case of location == 'centroid' the dimension values must
193        be a list of a Numerical array of length N, N being the number
194        of elements. Otherwise it must be of dimension Nx3
195
196        The returned values with be a list the length of indexes
197        (N if indexes = None).  Each value will be a list of the three
198        vertex values for this quantity.
199
200        Indexes is the set of element ids that the operation applies to.
201
202        """
203        from Numeric import take
204
205        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
206            msg = 'Invalid location: %s' %location
207            raise msg
208
209        import types, Numeric
210        assert type(indexes) in [types.ListType, types.NoneType,
211                                 Numeric.ArrayType],\
212                                 'Indices must be a list or None'
213
214        if location == 'centroids':
215            if (indexes ==  None):
216                indexes = range(len(self))
217            return take(self.centroid_values,indexes)
218        elif location == 'edges':
219            if (indexes ==  None):
220                indexes = range(len(self))
221            return take(self.edge_values,indexes)
222        elif location == 'unique vertices':
223            if (indexes ==  None):
224                indexes=range(self.domain.coordinates.shape[0])
225            vert_values = []
226            #Go through list of unique vertices
227            for unique_vert_id in indexes:
228                triangles = self.domain.vertexlist[unique_vert_id]
229
230                #In case there are unused points
231                if triangles is None:
232                    msg = 'Unique vertex not associated with triangles'
233                    raise msg
234
235                # Go through all triangle, vertex pairs
236                # Average the values
237                sum = 0
238                for triangle_id, vertex_id in triangles:
239                    sum += self.vertex_values[triangle_id, vertex_id]
240                vert_values.append(sum/len(triangles))
241            return Numeric.array(vert_values)
242        else:
243            if (indexes ==  None):
244                indexes = range(len(self))
245            return take(self.vertex_values,indexes)
246
247
248    def set_function_values(self, f, location='vertices', indexes = None):
249        """Set values for quantity using specified function
250
251        f: x, y -> z Function where x, y and z are arrays
252        location: Where values are to be stored.
253                  Permissible options are: vertices, centroid
254                  Default is "vertices"
255        """
[1626]256
257        #FIXME: Should check that function returns something sensible and
258        #raise a meaningfol exception if it returns None for example
259       
[1290]260        from Numeric import take
261
262        if (indexes ==  None):
263            indexes = range(len(self))
264            is_subset = False
265        else:
266            is_subset = True
267        if location == 'centroids':
268            P = take(self.domain.centroid_coordinates,indexes)
269            if is_subset:
270                self.set_values(f(P[:,0], P[:,1]), location, indexes = indexes)
271            else:
272                self.set_values(f(P[:,0], P[:,1]), location)
273        elif location == 'vertices':
274            P = self.domain.vertex_coordinates
275            if is_subset:
276                #Brute force
277                for e in indexes:
278                    for i in range(3):
279                        self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1])
280            else:
281                for i in range(3):
282                    self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
283        else:
284            raise 'Not implemented: %s' %location
285
286
287    def set_array_values(self, values, location='vertices', indexes = None):
288        """Set values for quantity
289
290        values: Numeric array
291        location: Where values are to be stored.
292        Permissible options are: vertices, edges, centroid, unique vertices
293        Default is "vertices"
294
295        indexes - if this action is carried out on a subset of
296        elements or unique vertices
297        The element/unique vertex indexes are specified here.
298
299        In case of location == 'centroid' the dimension values must
300        be a list of a Numerical array of length N, N being the number
301        of elements.
302
303        Otherwise it must be of dimension Nx3
304
305        The values will be stored in elements following their
306        internal ordering.
307
308        If selected location is vertices, values for centroid and edges
309        will be assigned interpolated values.
310        In any other case, only values for the specified locations
311        will be assigned and the others will be left undefined.
312        """
313
314        from Numeric import array, Float, Int, allclose
315
316        values = array(values).astype(Float)
317
318        if (indexes <> None):
319            indexes = array(indexes).astype(Int)
320            msg = 'Number of values must match number of indexes'
321            assert values.shape[0] == indexes.shape[0], msg
322
323        N = self.centroid_values.shape[0]
324
325        if location == 'centroids':
326            assert len(values.shape) == 1, 'Values array must be 1d'
327
328            if indexes == None:
329                msg = 'Number of values must match number of elements'
330                assert values.shape[0] == N, msg
331
332                self.centroid_values = values
333            else:
334                msg = 'Number of values must match number of indexes'
335                assert values.shape[0] == indexes.shape[0], msg
336
337                #Brute force
338                for i in range(len(indexes)):
339                    self.centroid_values[indexes[i]] = values[i]
340
341        elif location == 'edges':
342            assert len(values.shape) == 2, 'Values array must be 2d'
343
344            msg = 'Number of values must match number of elements'
345            assert values.shape[0] == N, msg
346
347            msg = 'Array must be N x 3'
348            assert values.shape[1] == 3, msg
349
350            self.edge_values = values
351
352        elif location == 'unique vertices':
353            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
354                   'Values array must be 1d'
355
356            self.set_vertex_values(values.flat, indexes = indexes)
357        else:
358            if len(values.shape) == 1:
359                self.set_vertex_values(values, indexes = indexes)
360                #if indexes == None:
361                    #Values are being specified once for each unique vertex
362                #    msg = 'Number of values must match number of vertices'
363                #    assert values.shape[0] == self.domain.coordinates.shape[0], msg
364                 #   self.set_vertex_values(values)
365                #else:
366                #    for element_index, value in map(None, indexes, values):
367                #        self.vertex_values[element_index, :] = value
368
369            elif len(values.shape) == 2:
370                #Vertex values are given as a triplet for each triangle
371
372                msg = 'Array must be N x 3'
373                assert values.shape[1] == 3, msg
374
375                if indexes == None:
376                    self.vertex_values = values
377                else:
378                    for element_index, value in map(None, indexes, values):
379                        self.vertex_values[element_index] = value
380            else:
381                msg = 'Values array must be 1d or 2d'
382                raise msg
383
384
385    # FIXME have a get_vertex_values as well, so the 'stage' quantity can be
386    # set, based on the elevation
387    def set_vertex_values(self, A, indexes = None):
388        """Set vertex values for all unique vertices based on input array A
389        which has one entry per unique vertex, i.e.
390        one value for each row in array self.domain.coordinates or
391        one value for each row in vertexlist.
392
393        indexes is the list of vertex_id's that will be set.
394
395        Note: Functions not allowed
396        """
397
398        from Numeric import array, Float
399
400        #Assert that A can be converted to a Numeric array of appropriate dim
401        A = array(A, Float)
402
403        #print 'SHAPE A', A.shape
404        assert len(A.shape) == 1
405
406        if indexes == None:
407            assert A.shape[0] == self.domain.coordinates.shape[0]
408            vertex_list = range(A.shape[0])
409        else:
410            assert A.shape[0] == len(indexes)
411            vertex_list = indexes
412        #Go through list of unique vertices
413        for i_index,unique_vert_id in enumerate(vertex_list):
414            triangles = self.domain.vertexlist[unique_vert_id]
415
416            if triangles is None: continue #In case there are unused points
417
418            #Go through all triangle, vertex pairs
419            #touching vertex unique_vert_id and set corresponding vertex value
420            for triangle_id, vertex_id in triangles:
421                self.vertex_values[triangle_id, vertex_id] = A[i_index]
422
423        #Intialise centroid and edge_values
424        self.interpolate()
425
426    def smooth_vertex_values(self, value_array='field_values',
427                             precision = None):
428        """ Smooths field_values or conserved_quantities data.
429        TODO: be able to smooth individual fields
430        NOTE:  This function does not have a test.
431        FIXME: NOT DONE - do we need it?
432        FIXME: this function isn't called by anything.
433               Maybe it should be removed..-DSG
434        """
435
436        from Numeric import concatenate, zeros, Float, Int, array, reshape
437
438
439        A,V = self.get_vertex_values(xy=False,
440                                     value_array=value_array,
441                                     smooth = True,
442                                     precision = precision)
443
444        #Set some field values
445        for volume in self:
446            for i,v in enumerate(volume.vertices):
447                if value_array == 'field_values':
448                    volume.set_field_values('vertex', i, A[v,:])
449                elif value_array == 'conserved_quantities':
450                    volume.set_conserved_quantities('vertex', i, A[v,:])
451
452        if value_array == 'field_values':
453            self.precompute()
454        elif value_array == 'conserved_quantities':
455            Volume.interpolate_conserved_quantities()
456
457
458    #Method for outputting model results
459    #FIXME: Split up into geometric and numeric stuff.
460    #FIXME: Geometric (X,Y,V) should live in mesh.py
461    #FIXME: STill remember to move XY to mesh
462    def get_vertex_values(self,
463                          xy=True,
464                          smooth = None,
465                          precision = None,
466                          reduction = None):
467        """Return vertex values like an OBJ format
468
469        The vertex values are returned as one sequence in the 1D float array A.
470        If requested the coordinates will be returned in 1D arrays X and Y.
471
472        The connectivity is represented as an integer array, V, of dimension
473        M x 3, where M is the number of volumes. Each row has three indices
474        into the X, Y, A arrays defining the triangle.
475
476        if smooth is True, vertex values corresponding to one common
477        coordinate set will be smoothed according to the given
478        reduction operator. In this case vertex coordinates will be
479        de-duplicated.
480
481        If no smoothings is required, vertex coordinates and values will
482        be aggregated as a concatenation of values at
483        vertices 0, vertices 1 and vertices 2
484
485
486        Calling convention
487        if xy is True:
488           X,Y,A,V = get_vertex_values
489        else:
490           A,V = get_vertex_values
491
492        """
493
494        from Numeric import concatenate, zeros, Float, Int, array, reshape
495
496
497        if smooth is None:
498            smooth = self.domain.smooth
499
500        if precision is None:
501            precision = Float
502
503        if reduction is None:
504            reduction = self.domain.reduction
505
506        #Create connectivity
507
508        if smooth == True:
509
510            V = self.domain.get_vertices()
511            N = len(self.domain.vertexlist)
512            A = zeros(N, precision)
513
514            #Smoothing loop
515            for k in range(N):
516                L = self.domain.vertexlist[k]
517
518                #Go through all triangle, vertex pairs
519                #contributing to vertex k and register vertex value
520
521                if L is None: continue #In case there are unused points
522
523                contributions = []
524                for volume_id, vertex_id in L:
525                    v = self.vertex_values[volume_id, vertex_id]
526                    contributions.append(v)
527
528                A[k] = reduction(contributions)
529
530
531            if xy is True:
532                X = self.domain.coordinates[:,0].astype(precision)
533                Y = self.domain.coordinates[:,1].astype(precision)
534
535                return X, Y, A, V
536            else:
537                return A, V
538        else:
539            #Don't smooth
[1632]540            #obj machinery moved to general_mesh
[1290]541
542            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
[1393]543            # These vert_id's will relate to the verts created below
[1632]544            #m = len(self.domain)  #Number of volumes
545            #M = 3*m        #Total number of unique vertices
546            #V = reshape(array(range(M)).astype(Int), (m,3))
547           
548            V = self.domain.get_triangles(obj=True)
549            #FIXME use get_vertices, when ready
[1290]550
551            A = self.vertex_values.flat
552
553            #Do vertex coordinates
554            if xy is True:
555                C = self.domain.get_vertex_coordinates()
556
557                X = C[:,0:6:2].copy()
558                Y = C[:,1:6:2].copy()
559
560                return X.flat, Y.flat, A, V
561            else:
562                return A, V
563
564
565    def extrapolate_first_order(self):
566        """Extrapolate conserved quantities from centroid to
567        vertices for each volume using
568        first order scheme.
569        """
570
571        qc = self.centroid_values
572        qv = self.vertex_values
573
574        for i in range(3):
575            qv[:,i] = qc
576
577
578    def get_integral(self):
579        """Compute the integral of quantity across entire domain
580        """
581        integral = 0
582        for k in range(self.domain.number_of_elements):
583            area = self.domain.areas[k]
584            qc = self.centroid_values[k]
585            integral += qc*area
586
587        return integral
588
589
590class Conserved_quantity(Quantity):
591    """Class conserved quantity adds to Quantity:
592
593    boundary values, storage and method for updating, and
594    methods for (second order) extrapolation from centroid to vertices inluding
595    gradients and limiters
596    """
597
598    def __init__(self, domain, vertex_values=None):
599        Quantity.__init__(self, domain, vertex_values)
600
601        from Numeric import zeros, Float
602
603        #Allocate space for boundary values
604        L = len(domain.boundary)
605        self.boundary_values = zeros(L, Float)
606
607        #Allocate space for updates of conserved quantities by
608        #flux calculations and forcing functions
609
610        N = domain.number_of_elements
611        self.explicit_update = zeros(N, Float )
612        self.semi_implicit_update = zeros(N, Float )
613
614
615    def update(self, timestep):
616        #Call correct module function
617        #(either from this module or C-extension)
618        return update(self, timestep)
619
620
621    def compute_gradients(self):
622        #Call correct module function
623        #(either from this module or C-extension)
624        return compute_gradients(self)
625
626
627    def limit(self):
628        #Call correct module function
629        #(either from this module or C-extension)
630        limit(self)
631
632
633    def extrapolate_second_order(self):
634        #Call correct module function
635        #(either from this module or C-extension)
636        extrapolate_second_order(self)
637
638
639def update(quantity, timestep):
640    """Update centroid values based on values stored in
641    explicit_update and semi_implicit_update as well as given timestep
642
643    Function implementing forcing terms must take on argument
644    which is the domain and they must update either explicit
645    or implicit updates, e,g,:
646
647    def gravity(domain):
648        ....
649        domain.quantities['xmomentum'].explicit_update = ...
650        domain.quantities['ymomentum'].explicit_update = ...
651
652
653
654    Explicit terms must have the form
655
656        G(q, t)
657
658    and explicit scheme is
659
660       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
661
662
663    Semi implicit forcing terms are assumed to have the form
664
665       G(q, t) = H(q, t) q
666
667    and the semi implicit scheme will then be
668
669      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
670
671
672    """
673
674    from Numeric import sum, equal, ones, Float
675
676    N = quantity.centroid_values.shape[0]
677
678
679    #Divide H by conserved quantity to obtain G (see docstring above)
680
681
682    for k in range(N):
683        x = quantity.centroid_values[k]
684        if x == 0.0:
685            #FIXME: Is this right
686            quantity.semi_implicit_update[k] = 0.0
687        else:
688            quantity.semi_implicit_update[k] /= x
689
690    #Explicit updates
691    quantity.centroid_values += timestep*quantity.explicit_update
692
693    #Semi implicit updates
694    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
695
696    if sum(equal(denominator, 0.0)) > 0.0:
697        msg = 'Zero division in semi implicit update. Call Stephen :-)'
698        raise msg
699    else:
700        #Update conserved_quantities from semi implicit updates
701        quantity.centroid_values /= denominator
702
703
704def interpolate_from_vertices_to_edges(quantity):
705    """Compute edge values from vertex values using linear interpolation
706    """
707
708    for k in range(quantity.vertex_values.shape[0]):
709        q0 = quantity.vertex_values[k, 0]
710        q1 = quantity.vertex_values[k, 1]
711        q2 = quantity.vertex_values[k, 2]
712
713        quantity.edge_values[k, 0] = 0.5*(q1+q2)
714        quantity.edge_values[k, 1] = 0.5*(q0+q2)
715        quantity.edge_values[k, 2] = 0.5*(q0+q1)
716
717
718
719def extrapolate_second_order(quantity):
720    """Extrapolate conserved quantities from centroid to
721    vertices for each volume using
722    second order scheme.
723    """
724
725    a, b = quantity.compute_gradients()
726
727    X = quantity.domain.get_vertex_coordinates()
728    qc = quantity.centroid_values
729    qv = quantity.vertex_values
730
731    #Check each triangle
732    for k in range(quantity.domain.number_of_elements):
733        #Centroid coordinates
734        x, y = quantity.domain.centroid_coordinates[k]
735
736        #vertex coordinates
737        x0, y0, x1, y1, x2, y2 = X[k,:]
738
739        #Extrapolate
740        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
741        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
742        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
743
744
745def compute_gradients(quantity):
746    """Compute gradients of triangle surfaces defined by centroids of
747    neighbouring volumes.
748    If one edge is on the boundary, use own centroid as neighbour centroid.
749    If two or more are on the boundary, fall back to first order scheme.
750    """
751
752    from Numeric import zeros, Float
753    from util import gradient
754
755    centroid_coordinates = quantity.domain.centroid_coordinates
756    surrogate_neighbours = quantity.domain.surrogate_neighbours
757    centroid_values = quantity.centroid_values
758    number_of_boundaries = quantity.domain.number_of_boundaries
759
760    N = centroid_values.shape[0]
761
762    a = zeros(N, Float)
763    b = zeros(N, Float)
764
765    for k in range(N):
766        if number_of_boundaries[k] < 2:
767            #Two or three true neighbours
768
769            #Get indices of neighbours (or self when used as surrogate)
770            k0, k1, k2 = surrogate_neighbours[k,:]
771
772            #Get data
773            q0 = centroid_values[k0]
774            q1 = centroid_values[k1]
775            q2 = centroid_values[k2]
776
777            x0, y0 = centroid_coordinates[k0] #V0 centroid
778            x1, y1 = centroid_coordinates[k1] #V1 centroid
779            x2, y2 = centroid_coordinates[k2] #V2 centroid
780
781            #Gradient
782            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
783
784        elif number_of_boundaries[k] == 2:
785            #One true neighbour
786
787            #Get index of the one neighbour
788            for k0 in surrogate_neighbours[k,:]:
789                if k0 != k: break
790            assert k0 != k
791
792            k1 = k  #self
793
794            #Get data
795            q0 = centroid_values[k0]
796            q1 = centroid_values[k1]
797
798            x0, y0 = centroid_coordinates[k0] #V0 centroid
799            x1, y1 = centroid_coordinates[k1] #V1 centroid
800
801            #Gradient
[1486]802            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)         
[1290]803        else:
804            #No true neighbours -
805            #Fall back to first order scheme
806            pass
807
808
809    return a, b
810
811
812
813def limit(quantity):
814    """Limit slopes for each volume to eliminate artificial variance
815    introduced by e.g. second order extrapolator
816
817    This is an unsophisticated limiter as it does not take into
818    account dependencies among quantities.
819
820    precondition:
821    vertex values are estimated from gradient
822    postcondition:
823    vertex values are updated
824    """
825
826    from Numeric import zeros, Float
827
828    N = quantity.domain.number_of_elements
829
830    beta_w = quantity.domain.beta_w
831
832    qc = quantity.centroid_values
833    qv = quantity.vertex_values
834
835    #Find min and max of this and neighbour's centroid values
836    qmax = zeros(qc.shape, Float)
837    qmin = zeros(qc.shape, Float)
838
839    for k in range(N):
840        qmax[k] = qmin[k] = qc[k]
841        for i in range(3):
842            n = quantity.domain.neighbours[k,i]
843            if n >= 0:
844                qn = qc[n] #Neighbour's centroid value
845
846                qmin[k] = min(qmin[k], qn)
847                qmax[k] = max(qmax[k], qn)
848
849
850    #Diffences between centroids and maxima/minima
851    dqmax = qmax - qc
852    dqmin = qmin - qc
853
854    #Deltas between vertex and centroid values
855    dq = zeros(qv.shape, Float)
856    for i in range(3):
857        dq[:,i] = qv[:,i] - qc
858
859    #Phi limiter
860    for k in range(N):
861
862        #Find the gradient limiter (phi) across vertices
863        phi = 1.0
864        for i in range(3):
865            r = 1.0
866            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
867            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
868
869            phi = min( min(r*beta_w, 1), phi )
870
871        #Then update using phi limiter
872        for i in range(3):
873            qv[k,i] = qc[k] + phi*dq[k,i]
874
875
876
877import compile
878if compile.can_use_C_extension('quantity_ext.c'):
879    #Replace python version with c implementations
880
881    from quantity_ext import limit, compute_gradients,\
882    extrapolate_second_order, interpolate_from_vertices_to_edges, update
Note: See TracBrowser for help on using the repository browser.