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

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

Did PEP8 pass, @brief.

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