source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py @ 4728

Last change on this file since 4728 was 4728, checked in by ole, 17 years ago

Cosmetics

File size: 52.7 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
17from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
18     argmax, argmin, allclose, take, reshape
19
20from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
21from anuga.utilities.polygon import inside_polygon
22
23from anuga.geospatial_data.geospatial_data import Geospatial_data
24from anuga.fit_interpolate.fit import fit_to_mesh
25from anuga.config import points_file_block_line_size as default_block_line_size
26
27class Quantity:
28
29    def __init__(self, domain, vertex_values=None):
30
31        from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
32
33        msg = 'First argument in Quantity.__init__ '
34        msg += 'must be of class Mesh (or a subclass thereof)'
35        assert isinstance(domain, Mesh), msg
36
37        if vertex_values is None:
38            N = len(domain) # number_of_elements
39            self.vertex_values = zeros((N, 3), Float)
40        else:
41            self.vertex_values = array(vertex_values).astype(Float)
42
43            N, V = self.vertex_values.shape
44            assert V == 3,\
45                   'Three vertex values per element must be specified'
46
47
48            msg = 'Number of vertex values (%d) must be consistent with'\
49                  %N
50            msg += 'number of elements in specified domain (%d).'\
51                   %len(domain)
52
53            assert N == len(domain), msg
54
55        self.domain = domain
56
57        #Allocate space for other quantities
58        self.centroid_values = zeros(N, Float)
59        self.edge_values = zeros((N, 3), Float)
60
61        #Intialise centroid and edge_values
62        self.interpolate()
63
64
65
66    #Methods for operator overloading
67    def __len__(self):
68        return self.centroid_values.shape[0]
69
70
71    def __neg__(self):
72        """Negate all values in this quantity giving meaning to the
73        expression -Q where Q is an instance of class Quantity
74        """
75
76        Q = Quantity(self.domain)
77        Q.set_values(-self.vertex_values)
78        return Q
79
80
81    def __add__(self, other):
82        """Add to self anything that could populate a quantity
83
84        E.g other can be a constant, an array, a function, another quantity
85        (except for a filename or points, attributes (for now))
86        - see set_values for details
87        """
88
89        Q = Quantity(self.domain)
90        Q.set_values(other)
91
92        result = Quantity(self.domain)
93        result.set_values(self.vertex_values + Q.vertex_values)
94        return result
95
96    def __radd__(self, other):
97        """Handle cases like 7+Q, where Q is an instance of class Quantity
98        """
99        return self + other
100
101
102    def __sub__(self, other):
103        return self + -other  #Invoke __neg__
104
105    def __mul__(self, other):
106        """Multiply self with anything that could populate a quantity
107
108        E.g other can be a constant, an array, a function, another quantity
109        (except for a filename or points, attributes (for now))
110        - see set_values for details
111
112        Note that if two quantitites q1 and q2 are multiplied,
113        vertex values are multiplied entry by entry
114        while centroid and edge values are re-interpolated.
115        Hence they won't be the product of centroid or edge values
116        from q1 and q2.
117        """
118
119        Q = Quantity(self.domain)
120        Q.set_values(other)
121
122        result = Quantity(self.domain)
123        result.set_values(self.vertex_values * Q.vertex_values)
124        return result
125
126    def __rmul__(self, other):
127        """Handle cases like 3*Q, where Q is an instance of class Quantity
128        """
129        return self * other
130
131    def __pow__(self, other):
132        """Raise quantity to (numerical) power
133
134        As with __mul__ vertex values are processed entry by entry
135        while centroid and edge values are re-interpolated.
136
137        Example using __pow__:
138          Q = (Q1**2 + Q2**2)**0.5
139
140        """
141
142        result = Quantity(self.domain)
143        result.set_values(self.vertex_values**other)
144        return result
145
146
147
148    def interpolate(self):
149        """Compute interpolated values at edges and centroid
150        Pre-condition: vertex_values have been set
151        """
152        # FIXME (Ole): This function is not called
153        # in real model runs. Consider removing it.
154       
155        N = self.vertex_values.shape[0]
156        for i in range(N):
157            v0 = self.vertex_values[i, 0]
158            v1 = self.vertex_values[i, 1]
159            v2 = self.vertex_values[i, 2]
160
161            self.centroid_values[i] = (v0 + v1 + v2)/3
162
163        self.interpolate_from_vertices_to_edges()
164
165
166    def interpolate_from_vertices_to_edges(self):
167        #Call correct module function
168        #(either from this module or C-extension)
169        interpolate_from_vertices_to_edges(self)
170
171
172
173
174    #New leaner interface to setting values
175    def set_values(self,
176                   numeric = None,    # List, numeric array or constant
177                   quantity = None,   # Another quantity
178                   function = None,   # Callable object: f(x,y)
179                   geospatial_data = None, # Arbitrary dataset
180                   points = None, values = None, data_georef = None, #Input
181                   # for fit (obsoleted by use of geo_spatial object)
182                   filename = None, attribute_name = None, #Input from file
183                   alpha = None,
184                   location = 'vertices',
185                   polygon = None,
186                   indices = None,
187                   verbose = False,
188                   use_cache = False):
189
190        """Set values for quantity based on different sources.
191
192        numeric:
193          Compatible list, Numeric array (see below) or constant.
194          If callable it will treated as a function (see below)
195          If instance of another Quantity it will be treated as such.
196          If geo_spatial object it will be treated as such
197
198        quantity:
199          Another quantity (compatible quantity, e.g. obtained as a
200          linear combination of quantities)
201
202        function:
203          Any callable object that takes two 1d arrays x and y
204          each of length N and returns an array also of length N.
205          The function will be evaluated at points determined by
206          location and indices in the underlying mesh.
207
208        geospatial_data:
209          Arbitrary geo spatial dataset in the form of the class
210          Geospatial_data. Mesh points are populated using
211          fit_interpolate.fit fitting
212
213        points:
214          Nx2 array of data points for use with fit_interpolate.fit
215          If points are present, an N array of attribute
216          values corresponding to
217          each data point must be present.
218          (Obsoleted by geospatial_data)         
219
220        values:
221          If points is specified, values is an array of length N containing
222          attribute values for each point.
223          (Obsoleted by geospatial_data)         
224
225        data_georef:
226          If points is specified, geo_reference applies to each point.
227          (Obsoleted by geospatial_data)         
228
229        filename:
230          Name of a points file containing data points and attributes for
231          use with fit_interpolate.fit.
232
233        attribute_name:
234          If specified, any array matching that name
235          will be used. from file or geospatial_data.
236          Otherwise a default will be used.
237
238        alpha:
239          Smoothing parameter to be used with fit_interpolate.fit.
240          See module fit_interpolate.fit for further details about alpha.
241          Alpha will only be used with points, values or filename.
242          Otherwise it will be ignored.
243
244
245        location: Where values are to be stored.
246                  Permissible options are: vertices, edges, centroids
247                  Default is 'vertices'
248
249                  In case of location == 'centroids' the dimension values must
250                  be a list of a Numerical array of length N,
251                  N being the number of elements.
252                  Otherwise it must be of dimension Nx3
253
254
255                  The values will be stored in elements following their
256                  internal ordering.
257
258                  If location is not 'unique vertices' Indices is the
259                  set of element ids that the operation applies to.
260                  If location is 'unique vertices' Indices is the set
261                  of vertex ids that the operation applies to.
262
263                  If selected location is vertices, values for
264                  centroid and edges will be assigned interpolated
265                  values.  In any other case, only values for the
266                  specified locations will be assigned and the others
267                  will be left undefined.
268
269
270        polygon: Restrict update of quantity to locations that fall
271                 inside polygon. Polygon works by selecting indices
272                 and calling set_values recursively.
273                 Polygon mode has only been implemented for
274                 constant values so far.
275
276        indices: Restrict update of quantity to locations that are
277                 identified by indices (e.g. node ids if location
278                 is 'vertices')       
279       
280        verbose: True means that output to stdout is generated
281
282        use_cache: True means that caching of intermediate results is
283                   attempted for fit_interpolate.fit.
284
285
286
287
288        Exactly one of the arguments
289          numeric, quantity, function, points, filename
290        must be present.
291        """
292
293        from anuga.geospatial_data.geospatial_data import Geospatial_data
294        from types import FloatType, IntType, LongType, ListType, NoneType
295        from Numeric import ArrayType
296
297
298        # Treat special case: Polygon situation
299        # Location will be ignored and set to 'centroids'
300        # FIXME (Ole): This needs to be generalised and
301        # perhaps the notion of location and indices simplified
302       
303        if polygon is not None:
304            if indices is not None:
305                msg = 'Only one of polygon and indices can be specified'
306                raise Exception, msg
307
308            msg = 'With polygon selected, set_quantity must provide '
309            msg += 'the keyword numeric and it must (currently) be '
310            msg += 'a constant.'
311            if numeric is None:
312                raise Exception, msg           
313            else:
314                # Check that numeric is as constant
315                assert type(numeric) in [FloatType, IntType, LongType], msg
316
317
318            location = 'centroids'
319
320
321            points = self.domain.get_centroid_coordinates(absolute=True)
322            indices = inside_polygon(points, polygon)
323           
324            self.set_values_from_constant(numeric,
325                                          location, indices, verbose)
326           
327            self.extrapolate_first_order()
328            return
329       
330       
331
332
333
334
335        #General input checks
336        L = [numeric, quantity, function, geospatial_data, points, filename]
337        msg = 'Exactly one of the arguments '+\
338              'numeric, quantity, function, geospatial_data, points, '+\
339              'or filename must be present.'
340        assert L.count(None) == len(L)-1, msg
341
342
343        if location not in ['vertices', 'centroids', 'edges',
344                            'unique vertices']:
345            msg = 'Invalid location: %s' %location
346            raise Exception, msg
347
348
349        msg = 'Indices must be a list or None'
350        assert type(indices) in [ListType, NoneType, ArrayType], msg
351
352
353        if not(points is None and values is None and data_georef is None):
354            from warnings import warn
355
356            msg = 'Using points, values or data_georef with set_quantity '
357            msg += 'is obsolete. Please use a Geospatial_data object instead.'
358            warn(msg, DeprecationWarning)
359
360
361
362        #Determine which 'set_values_from_...' to use
363
364        if numeric is not None:
365            if type(numeric) in [FloatType, IntType, LongType]:
366                self.set_values_from_constant(numeric,
367                                              location, indices, verbose)
368            elif type(numeric) in [ArrayType, ListType]:
369                self.set_values_from_array(numeric,
370                                           location, indices, verbose)
371            elif callable(numeric):
372                self.set_values_from_function(numeric,
373                                              location, indices, verbose)
374            elif isinstance(numeric, Quantity):
375                self.set_values_from_quantity(numeric,
376                                              location, indices, verbose)
377            elif isinstance(numeric, Geospatial_data):
378                self.set_values_from_geospatial_data(numeric,
379                                                     alpha,
380                                                     location, indices,
381                                                     verbose = verbose,
382                                                     use_cache = use_cache)
383            else:
384                msg = 'Illegal type for argument numeric: %s' %str(numeric)
385                raise msg
386
387        elif quantity is not None:
388            self.set_values_from_quantity(quantity,
389                                          location, indices, verbose)
390        elif function is not None:
391            msg = 'Argument function must be callable'
392            assert callable(function), msg
393            self.set_values_from_function(function,
394                                          location, indices, verbose)
395        elif geospatial_data is not None:
396                self.set_values_from_geospatial_data(geospatial_data,
397                                                     alpha,
398                                                     location, indices,
399                                                     verbose = verbose,
400                                                     use_cache = use_cache)
401        elif points is not None:
402            print 'The usage of points in set_values will be deprecated.' +\
403                  'Please use the geospatial_data object.'
404
405            msg = 'When points are specified, associated values must also be.'
406            assert values is not None, msg
407            self.set_values_from_points(points, values, alpha,
408                                        location, indices,
409                                        data_georef = data_georef,
410                                        verbose = verbose,
411                                        use_cache = use_cache)
412        elif filename is not None:
413            if hasattr(self.domain, 'points_file_block_line_size'):
414                max_read_lines = self.domain.points_file_block_line_size
415            else:
416                max_read_lines = default_block_line_size
417            self.set_values_from_file(filename, attribute_name, alpha,
418                                      location, indices,
419                                      verbose = verbose,
420                                      max_read_lines=max_read_lines,
421                                      use_cache = use_cache)
422        else:
423            raise Exception, 'This can\'t happen :-)'
424
425
426
427        # Update all locations in triangles
428        if location == 'vertices' or location == 'unique vertices':
429            # Intialise centroid and edge_values
430            self.interpolate()
431
432        if location == 'centroids':
433            # Extrapolate 1st order - to capture notion of area being specified
434            self.extrapolate_first_order()
435
436
437
438    #Specific functions for setting values
439    def set_values_from_constant(self, X,
440                                 location, indices, verbose):
441        """Set quantity values from specified constant X
442        """
443
444        # FIXME (Ole): Somehow indices refer to centroids
445        # rather than vertices as default. See unit test
446        # test_set_vertex_values_using_general_interface_with_subset(self):
447       
448
449        if location == 'centroids':
450            if indices is None:
451                self.centroid_values[:] = X
452            else:
453                #Brute force
454                for i in indices:
455                    self.centroid_values[i] = X
456
457        elif location == 'edges':
458            if indices is None:
459                self.edge_values[:] = X
460            else:
461                #Brute force
462                for i in indices:
463                    self.edge_values[i] = X
464
465        elif location == 'unique vertices':
466            if indices is None:
467                self.edge_values[:] = X  #FIXME (Ole): Shouldn't this be vertex_values?
468            else:
469
470                #Go through list of unique vertices
471                for unique_vert_id in indices:
472
473                    triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
474                   
475                    #In case there are unused points
476                    if len(triangles) == 0:
477                        continue
478                   
479                    #Go through all triangle, vertex pairs
480                    #and set corresponding vertex value
481                    for triangle_id, vertex_id in triangles:
482                        self.vertex_values[triangle_id, vertex_id] = X
483
484                    #Intialise centroid and edge_values
485                    self.interpolate()
486        else:
487            if indices is None:
488                self.vertex_values[:] = X
489            else:
490                #Brute force
491                for i_vertex in indices:
492                    self.vertex_values[i_vertex] = X
493
494
495
496
497    def set_values_from_array(self, values,
498                              location='vertices',
499                              indices=None,
500                              verbose=False):
501        """Set values for quantity
502
503        values: Numeric array
504        location: Where values are to be stored.
505        Permissible options are: vertices, edges, centroid, unique vertices
506        Default is 'vertices'
507
508        indices - if this action is carried out on a subset of
509        elements or unique vertices
510        The element/unique vertex indices are specified here.
511
512        In case of location == 'centroid' the dimension values must
513        be a list of a Numerical array of length N, N being the number
514        of elements.
515
516        Otherwise it must be of dimension Nx3
517
518        The values will be stored in elements following their
519        internal ordering.
520
521        If selected location is vertices, values for centroid and edges
522        will be assigned interpolated values.
523        In any other case, only values for the specified locations
524        will be assigned and the others will be left undefined.
525        """
526
527        from Numeric import array, Float, Int, allclose
528
529        values = array(values).astype(Float)
530
531        if indices is not None:
532            indices = array(indices).astype(Int)
533            msg = 'Number of values must match number of indices'
534            assert values.shape[0] == indices.shape[0], msg
535
536        N = self.centroid_values.shape[0]
537
538        if location == 'centroids':
539            assert len(values.shape) == 1, 'Values array must be 1d'
540
541            if indices is None:
542                msg = 'Number of values must match number of elements'
543                assert values.shape[0] == N, msg
544
545                self.centroid_values = values
546            else:
547                msg = 'Number of values must match number of indices'
548                assert values.shape[0] == indices.shape[0], msg
549
550                #Brute force
551                for i in range(len(indices)):
552                    self.centroid_values[indices[i]] = values[i]
553
554        elif location == 'edges':
555            # FIXME (Ole): No mention of indices here. However, I don't
556            # think we ever need to set values at edges anyway
557            assert len(values.shape) == 2, 'Values array must be 2d'
558
559            msg = 'Number of values must match number of elements'
560            assert values.shape[0] == N, msg
561
562            msg = 'Array must be N x 3'
563            assert values.shape[1] == 3, msg
564
565            self.edge_values = values
566
567        elif location == 'unique vertices':
568            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
569                   'Values array must be 1d'
570
571            self.set_vertex_values(values.flat, indices=indices)
572           
573        else:
574            # Location vertices
575            if len(values.shape) == 1:
576                self.set_vertex_values(values, indices=indices)
577
578            elif len(values.shape) == 2:
579                #Vertex values are given as a triplet for each triangle
580
581                msg = 'Array must be N x 3'
582                assert values.shape[1] == 3, msg
583
584                if indices is None:
585                    self.vertex_values = values
586                else:
587                    for element_index, value in map(None, indices, values):
588                        self.vertex_values[element_index] = value
589            else:
590                msg = 'Values array must be 1d or 2d'
591                raise msg
592           
593
594    def set_values_from_quantity(self, q,
595                                 location, indices, verbose):
596        """Set quantity values from specified quantity instance q
597
598        Location is ignored - vertices will always be used here.
599        """
600
601
602        A = q.vertex_values
603
604        from Numeric import allclose
605        msg = 'Quantities are defined on different meshes. '+\
606              'This might be a case for implementing interpolation '+\
607              'between different meshes.'
608        assert allclose(A.shape, self.vertex_values.shape), msg
609
610        self.set_values(A, location='vertices',
611                        indices=indices,
612                        verbose=verbose)
613
614
615    def set_values_from_function(self, f,
616                                 location='vertices',
617                                 indices=None,
618                                 verbose=False):
619        """Set values for quantity using specified function
620
621        Input
622       
623        f: x, y -> z Function where x, y and z are arrays
624        location: Where values are to be stored.
625                  Permissible options are: vertices, centroid, edges,
626                  unique vertices
627                  Default is "vertices"
628        indices: 
629
630                 
631        """
632
633        #FIXME: Should check that function returns something sensible and
634        #raise a meaningfull exception if it returns None for example
635
636        #FIXME: Should supply absolute coordinates
637
638
639        # Compute the function values and call set_values again
640        if location == 'centroids':
641            if indices is None:
642                indices = range(len(self))
643               
644            V = take(self.domain.get_centroid_coordinates(), indices)
645            self.set_values(f(V[:,0], V[:,1]),
646                            location=location,
647                            indices=indices)
648           
649        elif location == 'vertices':
650
651            M = self.domain.number_of_triangles
652            V = self.domain.get_vertex_coordinates()
653
654            x = V[:,0]; y = V[:,1];                     
655            values = f(x, y)
656
657
658            # FIXME (Ole): This code should replace all the
659            # rest of this function and it would work, except
660            # one unit test in test_region fails.
661            # If that could be resolved this one will be
662            # more robust and simple.
663           
664            #values = reshape(values, (M,3))
665            #self.set_values(values,
666            #                location='vertices',
667            #                indices=indices)
668
669
670            # This should be removed
671            if is_scalar(values):
672                # Function returned a constant value
673                self.set_values_from_constant(values,
674                                              location, indices, verbose)
675                return
676
677            # This should be removed           
678            if indices is None:
679                for j in range(3):
680                    self.vertex_values[:,j] = values[j::3]                 
681            else:   
682                #Brute force
683                for i in indices:
684                    for j in range(3):
685                        self.vertex_values[i,j] = values[3*i+j]
686
687
688        else:
689            raise 'Not implemented: %s' %location
690
691
692
693    def set_values_from_geospatial_data(self, geospatial_data, alpha,
694                                        location, indices,
695                                        verbose = False,
696                                        use_cache = False):
697
698        #FIXME: Use this function for the time being. Later move code in here
699
700        points = geospatial_data.get_data_points(absolute = False)
701        values = geospatial_data.get_attributes()
702        data_georef = geospatial_data.get_geo_reference()
703
704
705
706        self.set_values_from_points(points, values, alpha,
707                                    location, indices,
708                                    data_georef = data_georef,
709                                    verbose = verbose,
710                                    use_cache = use_cache)
711
712
713
714    def set_values_from_points(self, points, values, alpha,
715                               location, indices,
716                               data_georef = None,
717                               verbose = False,
718                               use_cache = False):
719        """
720        Set quantity values from arbitray data points using
721        fit_interpolate.fit
722        """
723
724
725        from anuga.coordinate_transforms.geo_reference import Geo_reference
726
727
728        points = ensure_numeric(points, Float)
729        values = ensure_numeric(values, Float)
730
731        if location != 'vertices':
732            msg = 'set_values_from_points is only defined for '+\
733                  'location=\'vertices\''
734            raise ms
735
736        coordinates = self.domain.get_nodes()
737        triangles = self.domain.triangles      #FIXME
738
739
740        #Take care of georeferencing
741        if data_georef is None:
742            data_georef = Geo_reference()
743
744
745        mesh_georef = self.domain.geo_reference
746
747        #print mesh_georef
748        #print data_georef
749        #print points
750
751
752        # Call fit_interpolate.fit function
753        args = (coordinates, triangles, points, values)
754        kwargs = {'data_origin': data_georef.get_origin(),
755                  'mesh_origin': mesh_georef.get_origin(),
756                  'alpha': alpha,
757                  'verbose': verbose}
758
759        vertex_attributes = apply(fit_to_mesh,
760                                  args, kwargs)       
761
762        # Call underlying method using array values
763        self.set_values_from_array(vertex_attributes,
764                                   location, indices, verbose)
765
766    def set_values_from_file(self, filename, attribute_name, alpha,
767                             location, indices,
768                             verbose = False,
769                             use_cache = False,
770                             max_read_lines=None):
771        """Set quantity based on arbitrary points in a points file
772        using attribute_name selects name of attribute
773        present in file.
774        If attribute_name is not specified, use first available attribute
775        as defined in geospatial_data.
776        """
777
778        from types import StringType
779        msg = 'Filename must be a text string'
780        assert type(filename) == StringType, msg
781
782
783        if location != 'vertices':
784            msg = 'set_values_from_points is only defined for '+\
785                  'location=\'vertices\''
786            raise msg
787
788        coordinates = self.domain.get_nodes(absolute=True)
789        triangles = self.domain.triangles      #FIXME
790           
791        vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
792                                        alpha=alpha,
793                                        attribute_name=attribute_name,
794                                        use_cache=use_cache,
795                                        verbose=verbose,
796                                        max_read_lines=max_read_lines)
797                                           
798        # Call underlying method using array values
799        self.set_values_from_array(vertex_attributes,
800                                   location, indices, verbose)
801
802   
803    def get_extremum_index(self, mode=None, indices=None):
804        """Return index for maximum or minimum value of quantity (on centroids)
805
806        Optional arguments:
807            mode is either 'max'(default) or 'min'.
808            indices is the set of element ids that the operation applies to.
809
810        Usage:
811            i = get_extreme_index()
812
813        Notes:
814            We do not seek the extremum at vertices as each vertex can
815            have multiple values - one for each triangle sharing it.
816
817            If there are multiple cells with same maximum value, the
818            first cell encountered in the triangle array is returned.
819        """
820
821        V = self.get_values(location='centroids', indices=indices)
822
823        # Always return absolute indices
824        if mode is None or mode == 'max':
825            i = argmax(V)
826        elif mode == 'min':   
827            i = argmin(V)
828
829           
830        if indices is None:
831            return i
832        else:
833            return indices[i]
834
835
836    def get_maximum_index(self, indices=None):
837        """See get extreme index for details
838        """
839
840        return self.get_extremum_index(mode='max',
841                                       indices=indices)
842
843
844       
845    def get_maximum_value(self, indices=None):
846        """Return maximum value of quantity (on centroids)
847
848        Optional argument:
849            indices is the set of element ids that the operation applies to.
850
851        Usage:
852            v = get_maximum_value()
853
854        Note, we do not seek the maximum at vertices as each vertex can
855        have multiple values - one for each triangle sharing it           
856        """
857
858
859        i = self.get_maximum_index(indices)
860        V = self.get_values(location='centroids') #, indices=indices)
861       
862        return V[i]
863       
864
865    def get_maximum_location(self, indices=None):
866        """Return location of maximum value of quantity (on centroids)
867
868        Optional argument:
869            indices is the set of element ids that the operation applies to.
870
871        Usage:
872            x, y = get_maximum_location()
873
874
875        Notes:
876            We do not seek the maximum at vertices as each vertex can
877            have multiple values - one for each triangle sharing it.
878
879            If there are multiple cells with same maximum value, the
880            first cell encountered in the triangle array is returned.       
881        """
882
883        i = self.get_maximum_index(indices)
884        x, y = self.domain.get_centroid_coordinates()[i]
885
886        return x, y
887
888
889    def get_minimum_index(self, indices=None):
890        """See get extreme index for details
891        """       
892
893        return self.get_extremum_index(mode='min',
894                                       indices=indices)
895
896
897    def get_minimum_value(self, indices=None):
898        """Return minimum value of quantity (on centroids)
899
900        Optional argument:
901            indices is the set of element ids that the operation applies to.
902
903        Usage:
904            v = get_minimum_value()
905
906        See get_maximum_value for more details.   
907        """
908
909
910        i = self.get_minimum_index(indices)
911        V = self.get_values(location='centroids')
912       
913        return V[i]
914       
915
916    def get_minimum_location(self, indices=None):
917        """Return location of minimum value of quantity (on centroids)
918
919        Optional argument:
920            indices is the set of element ids that the operation applies to.
921
922        Usage:
923            x, y = get_minimum_location()
924
925
926        Notes:
927            We do not seek the maximum at vertices as each vertex can
928            have multiple values - one for each triangle sharing it.
929
930            If there are multiple cells with same maximum value, the
931            first cell encountered in the triangle array is returned.       
932        """
933
934        i = self.get_minimum_index(indices)
935        x, y = self.domain.get_centroid_coordinates()[i]
936
937        return x, y
938
939
940
941
942    def get_interpolated_values(self, interpolation_points):
943
944        # Interpolation object based on internal (discontinuous triangles)
945        x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
946                                                                smooth=False)
947        # FIXME: This concat should roll into get_vertex_values
948        vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
949                                         axis=1)
950
951        can_reuse = False
952        if hasattr(self, 'interpolation_object'):
953            # Reuse to save time
954            I = self.interpolation_object
955
956            if allclose(interpolation_points, I._point_coordinates):
957                can_reuse = True
958               
959
960        if can_reuse is True:
961            # Use absence of points to indicate reuse in I.interpolate
962            result = I.interpolate(vertex_values) 
963        else:   
964            from anuga.fit_interpolate.interpolate import Interpolate
965
966            # Create interpolation object with matrix
967            I = Interpolate(vertex_coordinates, triangles)
968            self.interpolation_object = I
969
970            # Call interpolate with points the first time
971            interpolation_points = ensure_numeric(interpolation_points, Float)
972            result = I.interpolate(vertex_values, interpolation_points)     
973
974        return result
975
976
977    def get_values(self, interpolation_points=None,
978                   location='vertices',
979                   indices = None):
980        """get values for quantity
981
982        return X, Compatible list, Numeric array (see below)
983        interpolation_points: List of x, y coordinates where value is
984        sought (using interpolation). If points are given, values of
985        location and indices are ignored
986       
987        location: Where values are to be stored.
988                  Permissible options are: vertices, edges, centroids
989                  and unique vertices. Default is 'vertices'
990
991
992        The returned values with be a list the length of indices
993        (N if indices = None).
994
995        In case of location == 'centroids' the dimension of returned
996        values will be a list or a Numerical array of length N, N being
997        the number of elements.
998       
999        In case of location == 'vertices' or 'edges' the dimension of
1000        returned values will be of dimension Nx3
1001
1002        In case of location == 'unique vertices' the average value at
1003        each vertex will be returned and the dimension of returned values
1004        will be a 1d array of length "number of vertices"
1005       
1006        Indices is the set of element ids that the operation applies to.
1007
1008        The values will be stored in elements following their
1009        internal ordering.
1010
1011        """
1012        from Numeric import take
1013
1014        if interpolation_points is not None:
1015            return self.get_interpolated_values(interpolation_points)
1016       
1017       
1018
1019        if location not in ['vertices', 'centroids', 'edges',
1020                            'unique vertices']:
1021            msg = 'Invalid location: %s' %location
1022            raise msg
1023
1024        import types, Numeric
1025        assert type(indices) in [types.ListType, types.NoneType,
1026                                 Numeric.ArrayType],\
1027                                 'Indices must be a list or None'
1028
1029        if location == 'centroids':
1030            if (indices ==  None):
1031                indices = range(len(self))
1032            return take(self.centroid_values,indices)
1033        elif location == 'edges':
1034            if (indices ==  None):
1035                indices = range(len(self))
1036            return take(self.edge_values,indices)
1037        elif location == 'unique vertices':
1038            if (indices ==  None):
1039                indices=range(self.domain.number_of_nodes)
1040            vert_values = []
1041            #Go through list of unique vertices
1042            for unique_vert_id in indices:
1043                triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1044                   
1045                #In case there are unused points
1046                if len(triangles) == 0:
1047                    msg = 'Unique vertex not associated with triangles'
1048                    raise msg
1049
1050                # Go through all triangle, vertex pairs
1051                # Average the values
1052               
1053                # FIXME (Ole): Should we merge this with get_vertex_values
1054                sum = 0
1055                for triangle_id, vertex_id in triangles:
1056                    sum += self.vertex_values[triangle_id, vertex_id]
1057                vert_values.append(sum/len(triangles))
1058            return Numeric.array(vert_values)
1059        else:
1060            if (indices is None):
1061                indices = range(len(self))
1062            return take(self.vertex_values, indices)
1063
1064
1065
1066    def set_vertex_values(self, A, indices = None):
1067        """Set vertex values for all unique vertices based on input array A
1068        which has one entry per unique vertex, i.e.
1069        one value for each row in array self.domain.nodes.
1070
1071        indices is the list of vertex_id's that will be set.
1072
1073        This function is used by set_values_from_array
1074        """
1075
1076        from Numeric import array, Float
1077
1078        #Assert that A can be converted to a Numeric array of appropriate dim
1079        A = array(A, Float)
1080
1081        #print 'SHAPE A', A.shape
1082        assert len(A.shape) == 1
1083
1084        if indices is None:
1085            assert A.shape[0] == self.domain.get_nodes().shape[0]
1086            vertex_list = range(A.shape[0])
1087        else:
1088            assert A.shape[0] == len(indices)
1089            vertex_list = indices
1090
1091        #Go through list of unique vertices
1092       
1093        for i_index, unique_vert_id in enumerate(vertex_list):
1094
1095
1096            triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1097                   
1098            #In case there are unused points
1099            if len(triangles) == 0: continue
1100
1101            #Go through all triangle, vertex pairs
1102            #touching vertex unique_vert_id and set corresponding vertex value
1103            for triangle_id, vertex_id in triangles:
1104                self.vertex_values[triangle_id, vertex_id] = A[i_index]
1105
1106        #Intialise centroid and edge_values
1107        self.interpolate()
1108
1109
1110    def smooth_vertex_values(self, value_array='field_values',
1111                             precision = None):
1112        """ Smooths field_values or conserved_quantities data.
1113        TODO: be able to smooth individual fields
1114        NOTE:  This function does not have a test.
1115        FIXME: NOT DONE - do we need it?
1116        FIXME: this function isn't called by anything.
1117               Maybe it should be removed..-DSG
1118        """
1119
1120        from Numeric import concatenate, zeros, Float, Int, array, reshape
1121
1122
1123        A,V = self.get_vertex_values(xy=False,
1124                                     value_array=value_array,
1125                                     smooth = True,
1126                                     precision = precision)
1127
1128        #Set some field values
1129        for volume in self:
1130            for i,v in enumerate(volume.vertices):
1131                if value_array == 'field_values':
1132                    volume.set_field_values('vertex', i, A[v,:])
1133                elif value_array == 'conserved_quantities':
1134                    volume.set_conserved_quantities('vertex', i, A[v,:])
1135
1136        if value_array == 'field_values':
1137            self.precompute()
1138        elif value_array == 'conserved_quantities':
1139            Volume.interpolate_conserved_quantities()
1140
1141
1142    # Methods for outputting model results
1143    def get_vertex_values(self,
1144                          xy=True,
1145                          smooth=None,
1146                          precision=None):
1147        """Return vertex values like an OBJ format i.e. one value per node.
1148
1149        The vertex values are returned as one sequence in the 1D float array A.
1150        If requested the coordinates will be returned in 1D arrays X and Y.
1151
1152        The connectivity is represented as an integer array, V, of dimension
1153        Mx3, where M is the number of triangles. Each row has three indices
1154        defining the triangle and they correspond to elements in the arrays
1155        X, Y and A.
1156
1157        if smooth is True, vertex values corresponding to one common
1158        coordinate set will be smoothed by taking the average of vertex values for each node.
1159        In this case vertex coordinates will be
1160        de-duplicated corresponding to the original nodes as obtained from
1161        the method general_mesh.get_nodes()
1162
1163        If no smoothings is required, vertex coordinates and values will
1164        be aggregated as a concatenation of values at
1165        vertices 0, vertices 1 and vertices 2. This corresponds to
1166        the node coordinates obtained from the method
1167        general_mesh.get_vertex_coordinates()
1168
1169
1170        Calling convention
1171        if xy is True:
1172           X,Y,A,V = get_vertex_values
1173        else:
1174           A,V = get_vertex_values
1175
1176        """
1177
1178        from Numeric import concatenate, zeros, Float, Int, array, reshape
1179
1180
1181        if smooth is None:
1182            # Take default from domain
1183            smooth = self.domain.smooth
1184
1185        if precision is None:
1186            precision = Float
1187           
1188
1189        if smooth is True:
1190            # Ensure continuous vertex values by averaging
1191            # values at each node
1192           
1193            V = self.domain.get_triangles()
1194            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1195            A = zeros(N, Float)
1196            points = self.domain.get_nodes()           
1197           
1198            if 1:
1199                # Fast C version
1200                average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1201                                      ensure_numeric(self.domain.number_of_triangles_per_node),
1202                                      ensure_numeric(self.vertex_values),
1203                                      A)
1204                A = A.astype(precision)
1205            else:   
1206
1207                # Slow Python version
1208               
1209                current_node = 0
1210                k = 0 # Track triangles touching on node
1211                total = 0.0
1212                for index in self.domain.vertex_value_indices:
1213                    if current_node == N:
1214                        msg = 'Current node exceeding number of nodes (%d) ' %(N)
1215                        raise msg
1216                   
1217
1218                   
1219                    k += 1
1220                   
1221                    volume_id = index / 3
1222                    vertex_id = index % 3
1223                 
1224                    #assert V[volume_id, vertex_id] == current_node
1225               
1226                    v = self.vertex_values[volume_id, vertex_id]
1227                    total += v
1228
1229                    #print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1230                    if self.domain.number_of_triangles_per_node[current_node] == k:
1231                        A[current_node] = total/k
1232               
1233                   
1234                        # Move on to next node
1235                        total = 0.0
1236                        k = 0
1237                        current_node += 1
1238
1239
1240
1241        else:
1242            # Allow discontinuous vertex values
1243            V = self.domain.get_disconnected_triangles()
1244            points = self.domain.get_vertex_coordinates()
1245            A = self.vertex_values.flat.astype(precision)
1246
1247
1248        # Return   
1249        if xy is True:
1250            X = points[:,0].astype(precision)
1251            Y = points[:,1].astype(precision)
1252           
1253            return X, Y, A, V
1254        else:
1255            return A, V           
1256
1257
1258
1259    def extrapolate_first_order(self):
1260        """Extrapolate conserved quantities from centroid to
1261        vertices for each volume using
1262        first order scheme.
1263        """
1264
1265        qc = self.centroid_values
1266        qv = self.vertex_values
1267
1268        for i in range(3):
1269            qv[:,i] = qc
1270
1271
1272    def get_integral(self):
1273        """Compute the integral of quantity across entire domain
1274        """
1275        integral = 0
1276        for k in range(len(self.domain)):
1277            area = self.domain.areas[k]
1278            qc = self.centroid_values[k]
1279            integral += qc*area
1280
1281        return integral
1282
1283
1284
1285
1286class Conserved_quantity(Quantity):
1287    """Class conserved quantity adds to Quantity:
1288
1289    boundary values, storage and method for updating, and
1290    methods for (second order) extrapolation from centroid to vertices inluding
1291    gradients and limiters
1292    """
1293
1294    def __init__(self, domain, vertex_values=None):
1295        Quantity.__init__(self, domain, vertex_values)
1296
1297        from Numeric import zeros, Float
1298
1299        #Allocate space for boundary values
1300        L = len(domain.boundary)
1301        self.boundary_values = zeros(L, Float)
1302
1303        #Allocate space for updates of conserved quantities by
1304        #flux calculations and forcing functions
1305
1306        N = len(domain) # number_of_triangles
1307        self.explicit_update = zeros(N, Float )
1308        self.semi_implicit_update = zeros(N, Float )
1309        self.centroid_backup_values = zeros(N, Float)       
1310
1311       
1312
1313
1314    def update(self, timestep):
1315        #Call correct module function
1316        #(either from this module or C-extension)
1317        return update(self, timestep)
1318
1319
1320    def compute_gradients(self):
1321        #Call correct module function
1322        #(either from this module or C-extension)
1323        return compute_gradients(self)
1324
1325    def limit(self):
1326        #Call correct module function
1327        #(either from this module or C-extension)
1328        limit(self)
1329
1330    def extrapolate_second_order(self):
1331        #Call correct module function
1332        #(either from this module or C-extension)
1333        extrapolate_second_order(self)
1334
1335    def backup_centroid_values(self):
1336        #Call correct module function
1337        #(either from this module or C-extension)
1338        backup_centroid_values(self)
1339
1340    def saxpy_centroid_values(self,a,b):
1341        #Call correct module function
1342        #(either from this module or C-extension)
1343        saxpy_centroid_values(self,a,b)
1344   
1345
1346def update(quantity, timestep):
1347    """Update centroid values based on values stored in
1348    explicit_update and semi_implicit_update as well as given timestep
1349
1350    Function implementing forcing terms must take on argument
1351    which is the domain and they must update either explicit
1352    or implicit updates, e,g,:
1353
1354    def gravity(domain):
1355        ....
1356        domain.quantities['xmomentum'].explicit_update = ...
1357        domain.quantities['ymomentum'].explicit_update = ...
1358
1359
1360
1361    Explicit terms must have the form
1362
1363        G(q, t)
1364
1365    and explicit scheme is
1366
1367       q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
1368
1369
1370    Semi implicit forcing terms are assumed to have the form
1371
1372       G(q, t) = H(q, t) q
1373
1374    and the semi implicit scheme will then be
1375
1376      q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
1377
1378
1379    """
1380
1381    from Numeric import sum, equal, ones, exp, Float
1382
1383    N = quantity.centroid_values.shape[0]
1384
1385
1386    # Divide H by conserved quantity to obtain G (see docstring above)
1387
1388
1389    for k in range(N):
1390        x = quantity.centroid_values[k]
1391        if x == 0.0:
1392            #FIXME: Is this right
1393            quantity.semi_implicit_update[k] = 0.0
1394        else:
1395            quantity.semi_implicit_update[k] /= x
1396
1397
1398    # Semi implicit updates
1399    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1400
1401    if sum(less(denominator, 1.0)) > 0.0:
1402        msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
1403        raise msg
1404
1405    if sum(equal(denominator, 0.0)) > 0.0:
1406        msg = 'Zero division in semi implicit update. Call Stephen :-)'
1407        raise msg
1408    else:
1409        #Update conserved_quantities from semi implicit updates
1410        quantity.centroid_values /= denominator
1411
1412#    quantity.centroid_values = exp(timestep*quantity.semi_implicit_update)*quantity.centroid_values
1413
1414    #Explicit updates
1415    quantity.centroid_values += timestep*quantity.explicit_update
1416
1417def interpolate_from_vertices_to_edges(quantity):
1418    """Compute edge values from vertex values using linear interpolation
1419    """
1420
1421    for k in range(quantity.vertex_values.shape[0]):
1422        q0 = quantity.vertex_values[k, 0]
1423        q1 = quantity.vertex_values[k, 1]
1424        q2 = quantity.vertex_values[k, 2]
1425
1426        quantity.edge_values[k, 0] = 0.5*(q1+q2)
1427        quantity.edge_values[k, 1] = 0.5*(q0+q2)
1428        quantity.edge_values[k, 2] = 0.5*(q0+q1)
1429
1430
1431
1432def backup_centroid_values(quantity):
1433    """Copy centroid values to backup array"""
1434
1435    qc = quantity.centroid_values
1436    qb = quantity.centroid_backup_values
1437
1438    #Check each triangle
1439    for k in range(len(quantity.domain)):
1440        qb[k] = qc[k]
1441
1442
1443def saxpy_centroid_values(quantity,a,b):
1444    """saxpy operation between centroid value and backup"""
1445
1446    qc = quantity.centroid_values
1447    qb = quantity.centroid_backup_values
1448
1449    #Check each triangle
1450    for k in range(len(quantity.domain)):
1451        qc[k] = a*qc[k]+b*qb[k]       
1452
1453
1454def compute_gradients(quantity):
1455    """Compute gradients of triangle surfaces defined by centroids of
1456    neighbouring volumes.
1457    If one edge is on the boundary, use own centroid as neighbour centroid.
1458    If two or more are on the boundary, fall back to first order scheme.
1459    """
1460
1461    from Numeric import zeros, Float
1462    from utilitites.numerical_tools import gradient
1463
1464    centroid_coordinates = quantity.domain.centroid_coordinates
1465    surrogate_neighbours = quantity.domain.surrogate_neighbours
1466    centroid_values = quantity.centroid_values
1467    number_of_boundaries = quantity.domain.number_of_boundaries
1468
1469    N = centroid_values.shape[0]
1470
1471    a = zeros(N, Float)
1472    b = zeros(N, Float)
1473
1474    for k in range(N):
1475        if number_of_boundaries[k] < 2:
1476            #Two or three true neighbours
1477
1478            #Get indices of neighbours (or self when used as surrogate)
1479            k0, k1, k2 = surrogate_neighbours[k,:]
1480
1481            #Get data
1482            q0 = centroid_values[k0]
1483            q1 = centroid_values[k1]
1484            q2 = centroid_values[k2]
1485
1486            x0, y0 = centroid_coordinates[k0] #V0 centroid
1487            x1, y1 = centroid_coordinates[k1] #V1 centroid
1488            x2, y2 = centroid_coordinates[k2] #V2 centroid
1489
1490            #Gradient
1491            a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1492
1493        elif number_of_boundaries[k] == 2:
1494            #One true neighbour
1495
1496            #Get index of the one neighbour
1497            for k0 in surrogate_neighbours[k,:]:
1498                if k0 != k: break
1499            assert k0 != k
1500
1501            k1 = k  #self
1502
1503            #Get data
1504            q0 = centroid_values[k0]
1505            q1 = centroid_values[k1]
1506
1507            x0, y0 = centroid_coordinates[k0] #V0 centroid
1508            x1, y1 = centroid_coordinates[k1] #V1 centroid
1509
1510            #Gradient
1511            a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1512        else:
1513            #No true neighbours -
1514            #Fall back to first order scheme
1515            pass
1516
1517
1518    return a, b
1519
1520
1521
1522def limit(quantity):
1523    """Limit slopes for each volume to eliminate artificial variance
1524    introduced by e.g. second order extrapolator
1525
1526    This is an unsophisticated limiter as it does not take into
1527    account dependencies among quantities.
1528
1529    precondition:
1530    vertex values are estimated from gradient
1531    postcondition:
1532    vertex values are updated
1533    """
1534
1535    from Numeric import zeros, Float
1536
1537    N = quantity.domain.number_of_nodes
1538
1539    beta_w = quantity.domain.beta_w
1540
1541    qc = quantity.centroid_values
1542    qv = quantity.vertex_values
1543
1544    #Find min and max of this and neighbour's centroid values
1545    qmax = zeros(qc.shape, Float)
1546    qmin = zeros(qc.shape, Float)
1547
1548    for k in range(N):
1549        qmax[k] = qc[k]
1550        qmin[k] = qc[k]
1551        for i in range(3):
1552            n = quantity.domain.neighbours[k,i]
1553            if n >= 0:
1554                qn = qc[n] #Neighbour's centroid value
1555
1556                qmin[k] = min(qmin[k], qn)
1557                qmax[k] = max(qmax[k], qn)
1558        qmax[k] = min(qmax[k], 2.0*qc[k])
1559        qmin[k] = max(qmin[k], 0.5*qc[k])
1560
1561
1562    #Diffences between centroids and maxima/minima
1563    dqmax = qmax - qc
1564    dqmin = qmin - qc
1565
1566    #Deltas between vertex and centroid values
1567    dq = zeros(qv.shape, Float)
1568    for i in range(3):
1569        dq[:,i] = qv[:,i] - qc
1570
1571    #Phi limiter
1572    for k in range(N):
1573
1574        #Find the gradient limiter (phi) across vertices
1575        phi = 1.0
1576        for i in range(3):
1577            r = 1.0
1578            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1579            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1580
1581            phi = min( min(r*beta_w, 1), phi )
1582
1583        #Then update using phi limiter
1584        for i in range(3):
1585            qv[k,i] = qc[k] + phi*dq[k,i]
1586
1587
1588
1589from anuga.utilities import compile
1590if compile.can_use_C_extension('quantity_ext.c'):
1591    #Replace python version with c implementations
1592
1593    from quantity_ext import average_vertex_values, backup_centroid_values, \
1594         saxpy_centroid_values
1595
1596    from quantity_ext import compute_gradients, limit,\
1597    extrapolate_second_order, interpolate_from_vertices_to_edges, update
Note: See TracBrowser for help on using the repository browser.