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

Last change on this file since 4712 was 4712, checked in by steve, 17 years ago

Working on 2nd order time stepping

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