source: inundation/ga/storm_surge/pyvolution-parallel/quantity.py @ 1454

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