source: branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py @ 7176

Last change on this file since 7176 was 7176, checked in by rwilson, 15 years ago

Back-merge from Numeric to numpy.

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