source: inundation/pyvolution/quantity.py @ 1928

Last change on this file since 1928 was 1928, checked in by ole, 18 years ago

Made create_quantity_from_expression a method

File size: 40.1 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
57
58    #Methods for operator overloading   
59    def __len__(self):
60        return self.centroid_values.shape[0]
61
62
63    def __neg__(self):
64        """Negate all values in this quantity giving meaning to the
65        expression -Q where Q is an instance of class Quantity
66        """
67
68        Q = Quantity(self.domain)
69        Q.set_values(-self.vertex_values)
70        return Q
71   
72
73    def __add__(self, other):
74        """Add to self anything that could populate a quantity
75       
76        E.g other can be a constant, an array, a function, another quantity
77        (except for a filename or points, attributes (for now))
78        - see set_values for details
79        """
80
81        Q = Quantity(self.domain)
82        Q.set_values(other)       
83
84        result = Quantity(self.domain)
85        result.set_values(self.vertex_values + Q.vertex_values)
86        return result
87
88    def __radd__(self, other):
89        """Handle cases like 7+Q, where Q is an instance of class Quantity
90        """
91        return self + other
92   
93   
94    def __sub__(self, other):   
95        return self + -other  #Invoke __neg__
96
97    def __mul__(self, other):
98        """Multiply self with anything that could populate a quantity
99       
100        E.g other can be a constant, an array, a function, another quantity
101        (except for a filename or points, attributes (for now))
102        - see set_values for details
103
104        Note that if two quantitites q1 and q2 are multiplied,
105        vertex values are multiplied entry by entry
106        while centroid and edge values are re-interpolated.
107        Hence they won't be the product of centroid or edge values
108        from q1 and q2.
109        """
110
111        Q = Quantity(self.domain)
112        Q.set_values(other)
113       
114        result = Quantity(self.domain)
115        result.set_values(self.vertex_values * Q.vertex_values)
116        return result
117       
118    def __rmul__(self, other):
119        """Handle cases like 3*Q, where Q is an instance of class Quantity
120        """
121        return self * other
122
123    def __pow__(self, other):
124        """Raise quantity to (numerical) power
125
126        As with __mul__ vertex values are processed entry by entry
127        while centroid and edge values are re-interpolated.
128
129        Example using __pow__:
130          Q = (Q1**2 + Q2**2)**0.5
131
132        """
133
134        result = Quantity(self.domain)
135        result.set_values(self.vertex_values**other)
136        return result
137       
138   
139
140    def interpolate(self):
141        """Compute interpolated values at edges and centroid
142        Pre-condition: vertex_values have been set
143        """
144
145        N = self.vertex_values.shape[0]
146        for i in range(N):
147            v0 = self.vertex_values[i, 0]
148            v1 = self.vertex_values[i, 1]
149            v2 = self.vertex_values[i, 2]
150
151            self.centroid_values[i] = (v0 + v1 + v2)/3
152
153        self.interpolate_from_vertices_to_edges()
154
155
156    def interpolate_from_vertices_to_edges(self):
157        #Call correct module function
158        #(either from this module or C-extension)
159        interpolate_from_vertices_to_edges(self)
160
161
162
163
164    #New leaner interface to setting values
165    def set_values(self,
166                   numeric = None,    # List, numeric array or constant
167                   quantity = None,   # Another quantity
168                   function = None,   # Callable object: f(x,y) 
169                   points = None, values = None, #Input for least squares
170                   filename = None, attribute_name = None, #Input from file
171                   alpha = None,
172                   location = 'vertices',                   
173                   indices = None,
174                   verbose = None,
175                   use_cache = False):
176       
177        """Set values for quantity based on different sources.
178       
179        numeric:
180          Compatible list, Numeric array (see below) or constant.
181          If callable it will treated as a function
182          If instance of another Quantity it will be treated as such.
183
184        quantity:
185          Another quantity (compatible quantity, e.g. obtained as a
186          linear combination of quantities)
187
188        function:
189          Any callable object that takes two 1d arrays x and y
190          each of length N and returns an array also of length N.
191          The function will be evaluated at points determined by
192          location and indices.
193
194        points:
195          Nx2 array of data points for use with least squares fit
196          If points are present, an N array of attribute
197          values corresponding to
198          each data point must be present.
199
200        filename:
201          Name of pts file containing data points and attributes for
202          use with least squares.
203          If attribute_name is specified, any array matching that name
204          will be used. Otherwise the first available one will be used.
205
206        alpha:
207          Smoothing parameter to be used with least squares fits.
208          See module least_squares for further details about alpha.
209          Alpha will only be used with points, values or filename.
210          Otherwise it will be ignored.
211               
212         
213        location: Where values are to be stored.
214                  Permissible options are: vertices, edges, centroids
215                  Default is 'vertices'
216
217                  In case of location == 'centroids' the dimension values must
218                  be a list of a Numerical array of length N,
219                  N being the number of elements.
220                  Otherwise it must be of dimension Nx3
221
222
223                  The values will be stored in elements following their
224                  internal ordering.
225
226                  If location is not 'unique vertices' Indices is the
227                  set of element ids that the operation applies to.
228                  If location is 'unique vertices' Indices is the set
229                  of vertex ids that the operation applies to.
230
231                  If selected location is vertices, values for
232                  centroid and edges will be assigned interpolated
233                  values.  In any other case, only values for the
234                  specified locations will be assigned and the others
235                  will be left undefined.
236
237        verbose: True means that output to stdout is generated
238
239        use_cache: True means that caching of intermediate results is
240                   attempted for least squares fit.
241                   
242                   
243       
244
245        Exactly one of the arguments
246          numeric, quantity, function, points, filename
247        must be present.
248        """
249
250        from types import FloatType, IntType, LongType, ListType, NoneType
251        from Numeric import ArrayType
252
253        #General input checks
254        L = [numeric, quantity, function, points, filename]       
255        msg = 'Exactly one of the arguments '+\
256              'numeric, quantity, function, points, or filename '+\
257              'must be present.'
258        assert L.count(None) == len(L)-1, msg
259
260
261        if location not in ['vertices', 'centroids', 'edges',
262                            'unique vertices']:
263            msg = 'Invalid location: %s' %location
264            raise msg
265
266
267        msg = 'Indices must be a list or None'
268        assert type(indices) in [ListType, NoneType, ArrayType], msg
269               
270
271       
272        #Determine which 'set_values_from_...' to use
273
274        if numeric is not None:
275            if type(numeric) in [FloatType, IntType, LongType]:
276                self.set_values_from_constant(numeric,
277                                              location, indices, verbose)
278            elif type(numeric) in [ArrayType, ListType]:
279                self.set_values_from_array(numeric,
280                                           location, indices, verbose)
281            elif callable(numeric):
282                self.set_values_from_function(numeric,
283                                              location, indices, verbose)
284            elif isinstance(numeric, Quantity):
285                self.set_values_from_quantity(numeric,
286                                              location, indices, verbose)
287            else:
288                msg = 'Illegal type for argument numeric: %s' %str(numeric)
289                raise msg
290           
291        elif quantity is not None:
292            self.set_values_from_quantity(quantity,
293                                          location, indices, verbose)
294        elif function is not None:
295            msg = 'Argument function must be callable'
296            assert callable(function), msg           
297            self.set_values_from_function(function,
298                                          location, indices, verbose)
299        elif points is not None:
300            msg = 'When points are specified, associated values must also be.'
301            assert values is not None, msg
302            self.set_values_from_points(points, values, alpha,
303                                        location, indices, verbose,
304                                        use_cache)
305        elif filename is not None:
306            self.set_values_from_file(filename, attribute_name, alpha,
307                                      location, indices, verbose,
308                                      use_cache)
309        else:
310            raise 'This can\'t happen :-)'
311
312
313        #Update all locations in triangles
314        if location == 'vertices' or location == 'unique vertices':
315            #Intialise centroid and edge_values
316            self.interpolate()
317
318        if location == 'centroids':
319            #Extrapolate 1st order - to capture notion of area being specified
320            self.extrapolate_first_order()
321
322               
323
324    #Specific functions for setting values
325    def set_values_from_constant(self, X,
326                                 location, indices, verbose):
327        """Set quantity values from specified constant X
328        """
329
330
331        if location == 'centroids':
332            if (indices ==  None):
333                self.centroid_values[:] = X
334            else:
335                #Brute force
336                for i in indices:
337                    self.centroid_values[i,:] = X
338
339        elif location == 'edges':
340            if (indices ==  None):
341                self.edge_values[:] = X
342            else:
343                #Brute force
344                for i in indices:
345                    self.edge_values[i,:] = X
346
347        elif location == 'unique vertices':
348            if (indices ==  None):
349                self.edge_values[:] = X
350            else:
351
352                #Go through list of unique vertices
353                for unique_vert_id in indices:
354                    triangles = self.domain.vertexlist[unique_vert_id]
355
356                    #In case there are unused points
357                    if triangles is None: continue
358
359                    #Go through all triangle, vertex pairs
360                    #and set corresponding vertex value
361                    for triangle_id, vertex_id in triangles:
362                        self.vertex_values[triangle_id, vertex_id] = X
363
364                    #Intialise centroid and edge_values
365                    self.interpolate()
366        else:
367            if (indices ==  None):
368                self.vertex_values[:] = X
369            else:
370                #Brute force
371                for i_vertex in indices:
372                    self.vertex_values[i_vertex,:] = X
373
374   
375
376   
377
378
379    def set_values_from_array(self, values,
380                              location, indices, verbose):
381        """Set values for quantity
382
383        values: Numeric array
384        location: Where values are to be stored.
385        Permissible options are: vertices, edges, centroid, unique vertices
386        Default is 'vertices'
387
388        indices - if this action is carried out on a subset of
389        elements or unique vertices
390        The element/unique vertex indices are specified here.
391
392        In case of location == 'centroid' the dimension values must
393        be a list of a Numerical array of length N, N being the number
394        of elements.
395
396        Otherwise it must be of dimension Nx3
397
398        The values will be stored in elements following their
399        internal ordering.
400
401        If selected location is vertices, values for centroid and edges
402        will be assigned interpolated values.
403        In any other case, only values for the specified locations
404        will be assigned and the others will be left undefined.
405        """
406
407        from Numeric import array, Float, Int, allclose
408
409        values = array(values).astype(Float)
410
411        if indices is not None:
412            indices = array(indices).astype(Int)
413            msg = 'Number of values must match number of indices'
414            assert values.shape[0] == indices.shape[0], msg
415
416        N = self.centroid_values.shape[0]
417
418        if location == 'centroids':
419            assert len(values.shape) == 1, 'Values array must be 1d'
420
421            if indices is None:
422                msg = 'Number of values must match number of elements'
423                assert values.shape[0] == N, msg
424
425                self.centroid_values = values
426            else:
427                msg = 'Number of values must match number of indices'
428                assert values.shape[0] == indices.shape[0], msg
429
430                #Brute force
431                for i in range(len(indices)):
432                    self.centroid_values[indices[i]] = values[i]
433
434        elif location == 'edges':
435            assert len(values.shape) == 2, 'Values array must be 2d'
436
437            msg = 'Number of values must match number of elements'
438            assert values.shape[0] == N, msg
439
440            msg = 'Array must be N x 3'
441            assert values.shape[1] == 3, msg
442
443            self.edge_values = values
444
445        elif location == 'unique vertices':
446            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
447                   'Values array must be 1d'
448
449            self.set_vertex_values(values.flat, indices = indices)
450        else:
451            if len(values.shape) == 1:
452                self.set_vertex_values(values, indices = indices)
453                #if indices == None:
454                    #Values are being specified once for each unique vertex
455                #    msg = 'Number of values must match number of vertices'
456                #    assert values.shape[0] == self.domain.coordinates.shape[0], msg
457                 #   self.set_vertex_values(values)
458                #else:
459                #    for element_index, value in map(None, indices, values):
460                #        self.vertex_values[element_index, :] = value
461
462            elif len(values.shape) == 2:
463                #Vertex values are given as a triplet for each triangle
464
465                msg = 'Array must be N x 3'
466                assert values.shape[1] == 3, msg
467
468                if indices == None:
469                    self.vertex_values = values
470                else:
471                    for element_index, value in map(None, indices, values):
472                        self.vertex_values[element_index] = value
473            else:
474                msg = 'Values array must be 1d or 2d'
475                raise msg
476
477    def set_values_from_quantity(self, q,
478                                 location, indices, verbose):
479        """Set quantity values from specified quantity instance q
480
481        Location is ignored
482        """
483       
484       
485        A = q.vertex_values
486
487        from Numeric import allclose
488        msg = 'Quantities are defined on different meshes. '+\
489              'This might be a case for implementing interpolation '+\
490              'between different meshes.'
491        assert allclose(A.shape, self.vertex_values.shape), msg
492
493        self.set_values(A, location='vertices',
494                        indices=indices,
495                        verbose=verbose)
496
497
498    def set_values_from_function(self, f,
499                                 location, indices, verbose):
500        """Set values for quantity using specified function
501
502        f: x, y -> z Function where x, y and z are arrays
503        location: Where values are to be stored.
504                  Permissible options are: vertices, centroid, edges,
505                  unique vertices
506                  Default is "vertices"
507        """
508
509        #FIXME: Should check that function returns something sensible and
510        #raise a meaningfull exception if it returns None for example
511
512        from Numeric import take
513
514        if (indices is None):
515            indices = range(len(self))
516            is_subset = False
517        else:
518            is_subset = True
519           
520        if location == 'centroids':
521            P = take(self.domain.centroid_coordinates, indices)
522            if is_subset:
523                self.set_values(f(P[:,0], P[:,1]),
524                                location = location,
525                                indices = indices)
526            else:
527                self.set_values(f(P[:,0], P[:,1]), location = location)
528        elif location == 'vertices':
529            P = self.domain.vertex_coordinates
530            if is_subset:
531                #Brute force
532                for e in indices:
533                    for i in range(3):
534                        self.vertex_values[e,i] = f(P[e,2*i], P[e,2*i+1])
535            else:
536                for i in range(3):
537                    self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
538        else:
539            raise 'Not implemented: %s' %location
540
541
542
543    def set_values_from_points(self, points, values, alpha,
544                               location, indices, verbose, use_cache):
545        """Set quantity values from arbitray data points using least squares
546        """
547
548        from Numeric import Float
549        from util import ensure_numeric
550        from least_squares import fit_to_mesh
551       
552        points = ensure_numeric(points, Float)
553        values = ensure_numeric(values, Float)
554
555        if location != 'vertices':
556            msg = 'set_values_from_points is only defined for'+\
557                  'location=\'vertices\''
558            raise msg
559
560        coordinates = self.domain.coordinates
561        triangles = self.domain.triangles
562
563        if use_cache is True:
564            try:
565                from caching import cache
566            except:
567                msg = 'Caching was requested, but caching module'+\
568                      'could not be imported'
569                raise msg
570           
571            args = (coordinates, triangles, points, values)
572            kwargs = {'alpha': alpha, 'verbose': verbose}
573            vertex_attributes = cache(fit_to_mesh,
574                                      args, kwargs,
575                                      verbose = verbose)
576        else:   
577            vertex_attributes = fit_to_mesh(coordinates,
578                                            triangles,
579                                            points,
580                                            values,
581                                            alpha = alpha,
582                                            verbose = verbose)
583           
584       
585        self.set_values_from_array(vertex_attributes,
586                                   location, indices, verbose)       
587       
588
589       
590       
591   
592    def set_values_from_file(self, filename, attribute_name, alpha,
593                             location, indices, verbose, use_cache):
594        """Set quantity based on arbitrary points in .pts file
595        using least_squares attribute_name selects name of attribute
596        present in file.
597        If not specified try to use whatever is available in file.
598        """
599
600
601        from types import StringType
602        msg = 'Filename must be a text string'
603        assert type(filename) == StringType, msg
604
605
606        #Read from (NetCDF) file
607        from load_mesh.loadASCII import import_points_file
608        points_dict = import_points_file(filename)
609        points = points_dict['pointlist']
610        attributes = points_dict['attributelist']
611       
612        if attribute_name is None:
613            names = attributes.keys()
614            attribute_name = names[0]
615
616        msg = 'Attribute_name must be a text string'
617        assert type(attribute_name) == StringType, msg           
618
619
620        if verbose:
621            print 'Using attribute %s from file %s' %(attribute_name, filename)
622            print 'Available attributes: %s' %(names)
623
624        try:
625            z = attributes[attribute_name]
626        except:
627            msg = 'Could not extract attribute %s from file %s'\
628                  %(attribute_name, filename)
629            raise msg
630
631           
632        #Call least squares method   
633        self.set_values_from_points(points, z, alpha,
634                                    location, indices, verbose, use_cache)
635
636
637
638    def get_values(self, location='vertices', indices = None):
639        """get values for quantity
640
641        return X, Compatible list, Numeric array (see below)
642        location: Where values are to be stored.
643                  Permissible options are: vertices, edges, centroid
644                  and unique vertices. Default is 'vertices'
645
646        In case of location == 'centroids' the dimension values must
647        be a list of a Numerical array of length N, N being the number
648        of elements. Otherwise it must be of dimension Nx3
649
650        The returned values with be a list the length of indices
651        (N if indices = None).  Each value will be a list of the three
652        vertex values for this quantity.
653
654        Indices is the set of element ids that the operation applies to.
655
656        """
657        from Numeric import take
658
659        if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
660            msg = 'Invalid location: %s' %location
661            raise msg
662
663        import types, Numeric
664        assert type(indices) in [types.ListType, types.NoneType,
665                                 Numeric.ArrayType],\
666                                 'Indices must be a list or None'
667
668        if location == 'centroids':
669            if (indices ==  None):
670                indices = range(len(self))
671            return take(self.centroid_values,indices)
672        elif location == 'edges':
673            if (indices ==  None):
674                indices = range(len(self))
675            return take(self.edge_values,indices)
676        elif location == 'unique vertices':
677            if (indices ==  None):
678                indices=range(self.domain.coordinates.shape[0])
679            vert_values = []
680            #Go through list of unique vertices
681            for unique_vert_id in indices:
682                triangles = self.domain.vertexlist[unique_vert_id]
683
684                #In case there are unused points
685                if triangles is None:
686                    msg = 'Unique vertex not associated with triangles'
687                    raise msg
688
689                # Go through all triangle, vertex pairs
690                # Average the values
691                sum = 0
692                for triangle_id, vertex_id in triangles:
693                    sum += self.vertex_values[triangle_id, vertex_id]
694                vert_values.append(sum/len(triangles))
695            return Numeric.array(vert_values)
696        else:
697            if (indices ==  None):
698                indices = range(len(self))
699            return take(self.vertex_values,indices)
700
701
702
703    def set_vertex_values(self, A, indices = None):
704        """Set vertex values for all unique vertices based on input array A
705        which has one entry per unique vertex, i.e.
706        one value for each row in array self.domain.coordinates or
707        one value for each row in vertexlist.
708
709        indices is the list of vertex_id's that will be set.
710
711        This function is used by set_values_from_array
712        """
713
714        from Numeric import array, Float
715
716        #Assert that A can be converted to a Numeric array of appropriate dim
717        A = array(A, Float)
718
719        #print 'SHAPE A', A.shape
720        assert len(A.shape) == 1
721
722        if indices == None:
723            assert A.shape[0] == self.domain.coordinates.shape[0]
724            vertex_list = range(A.shape[0])
725        else:
726            assert A.shape[0] == len(indices)
727            vertex_list = indices
728
729        #Go through list of unique vertices
730        for i_index, unique_vert_id in enumerate(vertex_list):
731            triangles = self.domain.vertexlist[unique_vert_id]
732
733            if triangles is None: continue #In case there are unused points
734
735            #Go through all triangle, vertex pairs
736            #touching vertex unique_vert_id and set corresponding vertex value
737            for triangle_id, vertex_id in triangles:
738                self.vertex_values[triangle_id, vertex_id] = A[i_index]
739
740        #Intialise centroid and edge_values
741        self.interpolate()
742
743
744    def smooth_vertex_values(self, value_array='field_values',
745                             precision = None):
746        """ Smooths field_values or conserved_quantities data.
747        TODO: be able to smooth individual fields
748        NOTE:  This function does not have a test.
749        FIXME: NOT DONE - do we need it?
750        FIXME: this function isn't called by anything.
751               Maybe it should be removed..-DSG
752        """
753
754        from Numeric import concatenate, zeros, Float, Int, array, reshape
755
756
757        A,V = self.get_vertex_values(xy=False,
758                                     value_array=value_array,
759                                     smooth = True,
760                                     precision = precision)
761
762        #Set some field values
763        for volume in self:
764            for i,v in enumerate(volume.vertices):
765                if value_array == 'field_values':
766                    volume.set_field_values('vertex', i, A[v,:])
767                elif value_array == 'conserved_quantities':
768                    volume.set_conserved_quantities('vertex', i, A[v,:])
769
770        if value_array == 'field_values':
771            self.precompute()
772        elif value_array == 'conserved_quantities':
773            Volume.interpolate_conserved_quantities()
774
775
776    #Method for outputting model results
777    #FIXME: Split up into geometric and numeric stuff.
778    #FIXME: Geometric (X,Y,V) should live in mesh.py
779    #FIXME: STill remember to move XY to mesh
780    def get_vertex_values(self,
781                          xy=True,
782                          smooth = None,
783                          precision = None,
784                          reduction = None):
785        """Return vertex values like an OBJ format
786
787        The vertex values are returned as one sequence in the 1D float array A.
788        If requested the coordinates will be returned in 1D arrays X and Y.
789
790        The connectivity is represented as an integer array, V, of dimension
791        M x 3, where M is the number of volumes. Each row has three indices
792        into the X, Y, A arrays defining the triangle.
793
794        if smooth is True, vertex values corresponding to one common
795        coordinate set will be smoothed according to the given
796        reduction operator. In this case vertex coordinates will be
797        de-duplicated.
798
799        If no smoothings is required, vertex coordinates and values will
800        be aggregated as a concatenation of values at
801        vertices 0, vertices 1 and vertices 2
802
803
804        Calling convention
805        if xy is True:
806           X,Y,A,V = get_vertex_values
807        else:
808           A,V = get_vertex_values
809
810        """
811
812        from Numeric import concatenate, zeros, Float, Int, array, reshape
813
814
815        if smooth is None:
816            smooth = self.domain.smooth
817
818        if precision is None:
819            precision = Float
820
821        if reduction is None:
822            reduction = self.domain.reduction
823
824        #Create connectivity
825
826        if smooth == True:
827
828            V = self.domain.get_vertices()
829            N = len(self.domain.vertexlist)
830            A = zeros(N, precision)
831
832            #Smoothing loop
833            for k in range(N):
834                L = self.domain.vertexlist[k]
835
836                #Go through all triangle, vertex pairs
837                #contributing to vertex k and register vertex value
838
839                if L is None: continue #In case there are unused points
840
841                contributions = []
842                for volume_id, vertex_id in L:
843                    v = self.vertex_values[volume_id, vertex_id]
844                    contributions.append(v)
845
846                A[k] = reduction(contributions)
847
848
849            if xy is True:
850                X = self.domain.coordinates[:,0].astype(precision)
851                Y = self.domain.coordinates[:,1].astype(precision)
852
853                return X, Y, A, V
854            else:
855                return A, V
856        else:
857            #Don't smooth
858            #obj machinery moved to general_mesh
859
860            # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]]
861            # These vert_id's will relate to the verts created below
862            #m = len(self.domain)  #Number of volumes
863            #M = 3*m        #Total number of unique vertices
864            #V = reshape(array(range(M)).astype(Int), (m,3))
865
866            V = self.domain.get_triangles(obj=True)
867            #FIXME use get_vertices, when ready
868
869            A = self.vertex_values.flat
870
871            #Do vertex coordinates
872            if xy is True:
873                C = self.domain.get_vertex_coordinates()
874
875                X = C[:,0:6:2].copy()
876                Y = C[:,1:6:2].copy()
877
878                return X.flat, Y.flat, A, V
879            else:
880                return A, V
881
882
883    def extrapolate_first_order(self):
884        """Extrapolate conserved quantities from centroid to
885        vertices for each volume using
886        first order scheme.
887        """
888
889        qc = self.centroid_values
890        qv = self.vertex_values
891
892        for i in range(3):
893            qv[:,i] = qc
894
895
896    def get_integral(self):
897        """Compute the integral of quantity across entire domain
898        """
899        integral = 0
900        for k in range(self.domain.number_of_elements):
901            area = self.domain.areas[k]
902            qc = self.centroid_values[k]
903            integral += qc*area
904
905        return integral
906
907
908   
909
910class Conserved_quantity(Quantity):
911    """Class conserved quantity adds to Quantity:
912
913    boundary values, storage and method for updating, and
914    methods for (second order) extrapolation from centroid to vertices inluding
915    gradients and limiters
916    """
917
918    def __init__(self, domain, vertex_values=None):
919        Quantity.__init__(self, domain, vertex_values)
920
921        from Numeric import zeros, Float
922
923        #Allocate space for boundary values
924        L = len(domain.boundary)
925        self.boundary_values = zeros(L, Float)
926
927        #Allocate space for updates of conserved quantities by
928        #flux calculations and forcing functions
929
930        N = domain.number_of_elements
931        self.explicit_update = zeros(N, Float )
932        self.semi_implicit_update = zeros(N, Float )
933
934
935    def update(self, timestep):
936        #Call correct module function
937        #(either from this module or C-extension)
938        return update(self, timestep)
939
940
941    def compute_gradients(self):
942        #Call correct module function
943        #(either from this module or C-extension)
944        return compute_gradients(self)
945
946
947    def limit(self):
948        #Call correct module function
949        #(either from this module or C-extension)
950        limit(self)
951
952
953    def extrapolate_second_order(self):
954        #Call correct module function
955        #(either from this module or C-extension)
956        extrapolate_second_order(self)
957
958
959def update(quantity, timestep):
960    """Update centroid values based on values stored in
961    explicit_update and semi_implicit_update as well as given timestep
962
963    Function implementing forcing terms must take on argument
964    which is the domain and they must update either explicit
965    or implicit updates, e,g,:
966
967    def gravity(domain):
968        ....
969        domain.quantities['xmomentum'].explicit_update = ...
970        domain.quantities['ymomentum'].explicit_update = ...
971
972
973
974    Explicit terms must have the form
975
976        G(q, t)
977
978    and explicit scheme is
979
980       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
981
982
983    Semi implicit forcing terms are assumed to have the form
984
985       G(q, t) = H(q, t) q
986
987    and the semi implicit scheme will then be
988
989      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
990
991
992    """
993
994    from Numeric import sum, equal, ones, Float
995
996    N = quantity.centroid_values.shape[0]
997
998
999    #Divide H by conserved quantity to obtain G (see docstring above)
1000
1001
1002    for k in range(N):
1003        x = quantity.centroid_values[k]
1004        if x == 0.0:
1005            #FIXME: Is this right
1006            quantity.semi_implicit_update[k] = 0.0
1007        else:
1008            quantity.semi_implicit_update[k] /= x
1009
1010    #Explicit updates
1011    quantity.centroid_values += timestep*quantity.explicit_update
1012
1013    #Semi implicit updates
1014    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1015
1016    if sum(equal(denominator, 0.0)) > 0.0:
1017        msg = 'Zero division in semi implicit update. Call Stephen :-)'
1018        raise msg
1019    else:
1020        #Update conserved_quantities from semi implicit updates
1021        quantity.centroid_values /= denominator
1022
1023
1024def interpolate_from_vertices_to_edges(quantity):
1025    """Compute edge values from vertex values using linear interpolation
1026    """
1027
1028    for k in range(quantity.vertex_values.shape[0]):
1029        q0 = quantity.vertex_values[k, 0]
1030        q1 = quantity.vertex_values[k, 1]
1031        q2 = quantity.vertex_values[k, 2]
1032
1033        quantity.edge_values[k, 0] = 0.5*(q1+q2)
1034        quantity.edge_values[k, 1] = 0.5*(q0+q2)
1035        quantity.edge_values[k, 2] = 0.5*(q0+q1)
1036
1037
1038
1039def extrapolate_second_order(quantity):
1040    """Extrapolate conserved quantities from centroid to
1041    vertices for each volume using
1042    second order scheme.
1043    """
1044
1045    a, b = quantity.compute_gradients()
1046
1047    X = quantity.domain.get_vertex_coordinates()
1048    qc = quantity.centroid_values
1049    qv = quantity.vertex_values
1050
1051    #Check each triangle
1052    for k in range(quantity.domain.number_of_elements):
1053        #Centroid coordinates
1054        x, y = quantity.domain.centroid_coordinates[k]
1055
1056        #vertex coordinates
1057        x0, y0, x1, y1, x2, y2 = X[k,:]
1058
1059        #Extrapolate
1060        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
1061        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
1062        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)
1063
1064
1065def compute_gradients(quantity):
1066    """Compute gradients of triangle surfaces defined by centroids of
1067    neighbouring volumes.
1068    If one edge is on the boundary, use own centroid as neighbour centroid.
1069    If two or more are on the boundary, fall back to first order scheme.
1070    """
1071
1072    from Numeric import zeros, Float
1073    from util import gradient
1074
1075    centroid_coordinates = quantity.domain.centroid_coordinates
1076    surrogate_neighbours = quantity.domain.surrogate_neighbours
1077    centroid_values = quantity.centroid_values
1078    number_of_boundaries = quantity.domain.number_of_boundaries
1079
1080    N = centroid_values.shape[0]
1081
1082    a = zeros(N, Float)
1083    b = zeros(N, Float)
1084
1085    for k in range(N):
1086        if number_of_boundaries[k] < 2:
1087            #Two or three true neighbours
1088
1089            #Get indices of neighbours (or self when used as surrogate)
1090            k0, k1, k2 = surrogate_neighbours[k,:]
1091
1092            #Get data
1093            q0 = centroid_values[k0]
1094            q1 = centroid_values[k1]
1095            q2 = centroid_values[k2]
1096
1097            x0, y0 = centroid_coordinates[k0] #V0 centroid
1098            x1, y1 = centroid_coordinates[k1] #V1 centroid
1099            x2, y2 = centroid_coordinates[k2] #V2 centroid
1100
1101            #Gradient
1102            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1103
1104        elif number_of_boundaries[k] == 2:
1105            #One true neighbour
1106
1107            #Get index of the one neighbour
1108            for k0 in surrogate_neighbours[k,:]:
1109                if k0 != k: break
1110            assert k0 != k
1111
1112            k1 = k  #self
1113
1114            #Get data
1115            q0 = centroid_values[k0]
1116            q1 = centroid_values[k1]
1117
1118            x0, y0 = centroid_coordinates[k0] #V0 centroid
1119            x1, y1 = centroid_coordinates[k1] #V1 centroid
1120
1121            #Gradient
1122            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1123        else:
1124            #No true neighbours -
1125            #Fall back to first order scheme
1126            pass
1127
1128
1129    return a, b
1130
1131
1132
1133def limit(quantity):
1134    """Limit slopes for each volume to eliminate artificial variance
1135    introduced by e.g. second order extrapolator
1136
1137    This is an unsophisticated limiter as it does not take into
1138    account dependencies among quantities.
1139
1140    precondition:
1141    vertex values are estimated from gradient
1142    postcondition:
1143    vertex values are updated
1144    """
1145
1146    from Numeric import zeros, Float
1147
1148    N = quantity.domain.number_of_elements
1149
1150    beta_w = quantity.domain.beta_w
1151
1152    qc = quantity.centroid_values
1153    qv = quantity.vertex_values
1154
1155    #Find min and max of this and neighbour's centroid values
1156    qmax = zeros(qc.shape, Float)
1157    qmin = zeros(qc.shape, Float)
1158
1159    for k in range(N):
1160        qmax[k] = qmin[k] = qc[k]
1161        for i in range(3):
1162            n = quantity.domain.neighbours[k,i]
1163            if n >= 0:
1164                qn = qc[n] #Neighbour's centroid value
1165
1166                qmin[k] = min(qmin[k], qn)
1167                qmax[k] = max(qmax[k], qn)
1168
1169
1170    #Diffences between centroids and maxima/minima
1171    dqmax = qmax - qc
1172    dqmin = qmin - qc
1173
1174    #Deltas between vertex and centroid values
1175    dq = zeros(qv.shape, Float)
1176    for i in range(3):
1177        dq[:,i] = qv[:,i] - qc
1178
1179    #Phi limiter
1180    for k in range(N):
1181
1182        #Find the gradient limiter (phi) across vertices
1183        phi = 1.0
1184        for i in range(3):
1185            r = 1.0
1186            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1187            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1188
1189            phi = min( min(r*beta_w, 1), phi )
1190
1191        #Then update using phi limiter
1192        for i in range(3):
1193            qv[k,i] = qc[k] + phi*dq[k,i]
1194
1195
1196
1197from utilities import compile
1198if compile.can_use_C_extension('quantity_ext.c'):
1199    #Replace python version with c implementations
1200
1201    from quantity_ext import limit, compute_gradients,\
1202    extrapolate_second_order, interpolate_from_vertices_to_edges, update
Note: See TracBrowser for help on using the repository browser.