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

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

numpy changes.

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