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

Last change on this file since 4768 was 4768, checked in by ole, 16 years ago

Backlog of comments and errormessages

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