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

Last change on this file since 5521 was 5521, checked in by ole, 14 years ago

Deprecated 'edges' as an option for location in set_values.
Validation tests pass as well as all unit tests except for two references to
this functionality.

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