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

Last change on this file since 4990 was 4990, checked in by steve, 16 years ago

Combining calculation of gradient, limiting edge by neighbour
and extrapolation in one procedure in quantity

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