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

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

Cleaned up according to styleguide etc

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