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

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

Added more keyword arguments and enabled caching of set_values_from_function.
There may be a problem if the function is a callable object, so this needs to be tested.

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