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

Last change on this file since 5520 was 5520, checked in by ole, 15 years ago

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