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

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

Comments and cleanup

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