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

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

Added better error messages and a new unit test that reveals
problem where fitting crashes by not being able to find
a triangle for a particular point within the mesh.

The test passes now, but if flag controlling mesh reuse
in quantity.py is switched on, it fails (on nautilus).

See ticket:314 for more info.

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