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

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

Added comment to quantity.py

File size: 58.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
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        if True:
878            # Use mesh as defined by domain
879            # This used to cause problems for caching due to quantities
880            # changing, but it now works using the appropriate Mesh object.
881            # This addressed ticket:242 and was made to work when bug
882            # in ticket:314 was fixed 18 March 2009.
883            vertex_attributes = fit_to_mesh(filename,
884                                            mesh=self.domain.mesh,
885                                            alpha=alpha,
886                                            attribute_name=attribute_name,
887                                            use_cache=use_cache,
888                                            verbose=verbose,
889                                            max_read_lines=max_read_lines)
890        else:
891            # This variant will cause Mesh object to be recreated
892            # in fit_to_mesh thus doubling up on the neighbour structure
893            # FIXME(Ole): This is now obsolete 19 Jan 2009 except for bug
894            # (ticket:314) which was fixed 18 March 2009.
895            nodes = self.domain.get_nodes(absolute=True)
896            triangles = self.domain.get_triangles()
897            vertex_attributes = fit_to_mesh(filename,
898                                            nodes, triangles,
899                                            mesh=None,
900                                            alpha=alpha,
901                                            attribute_name=attribute_name,
902                                            use_cache=use_cache,
903                                            verbose=verbose,
904                                            max_read_lines=max_read_lines)
905
906        # Call underlying method using array values
907        if verbose:
908            print 'Applying fitted data to domain'
909        self.set_values_from_array(vertex_attributes, location,
910                                   indices, use_cache=use_cache,
911                                   verbose=verbose)
912
913    ##
914    # @brief Get index for maximum or minimum value of quantity.
915    # @param mode Either 'max' or 'min'.
916    # @param indices Set of IDs of elements to work on.
917    def get_extremum_index(self, mode=None, indices=None):
918        """Return index for maximum or minimum value of quantity (on centroids)
919
920        Optional arguments:
921            mode is either 'max'(default) or 'min'.
922            indices is the set of element ids that the operation applies to.
923
924        Usage:
925            i = get_extreme_index()
926
927        Notes:
928            We do not seek the extremum at vertices as each vertex can
929            have multiple values - one for each triangle sharing it.
930
931            If there are multiple cells with same maximum value, the
932            first cell encountered in the triangle array is returned.
933        """
934
935        V = self.get_values(location='centroids', indices=indices)
936
937        # Always return absolute indices
938        if mode is None or mode == 'max':
939            i = num.argmax(V)
940        elif mode == 'min':
941            i = num.argmin(V)
942        else:
943            raise ValueError, 'Bad mode value, got: %s' % str(mode)
944
945        if indices is None:
946            return i
947        else:
948            return indices[i]
949
950    ##
951    # @brief Get index for maximum value of quantity.
952    # @param indices Set of IDs of elements to work on.
953    def get_maximum_index(self, indices=None):
954        """See get extreme index for details"""
955
956        return self.get_extremum_index(mode='max', indices=indices)
957
958    ##
959    # @brief Return maximum value of quantity (on centroids).
960    # @param indices Set of IDs of elements to work on.
961    def get_maximum_value(self, indices=None):
962        """Return maximum value of quantity (on centroids)
963
964        Optional argument:
965            indices is the set of element ids that the operation applies to.
966
967        Usage:
968            v = get_maximum_value()
969
970        Note, we do not seek the maximum at vertices as each vertex can
971        have multiple values - one for each triangle sharing it
972        """
973
974        i = self.get_maximum_index(indices)
975        V = self.get_values(location='centroids')      #, indices=indices)
976
977        return V[i]
978
979    ##
980    # @brief Get location of maximum value of quantity (on centroids).
981    # @param indices Set of IDs of elements to work on.
982    def get_maximum_location(self, indices=None):
983        """Return location of maximum value of quantity (on centroids)
984
985        Optional argument:
986            indices is the set of element ids that the operation applies to.
987
988        Usage:
989            x, y = get_maximum_location()
990
991        Notes:
992            We do not seek the maximum at vertices as each vertex can
993            have multiple values - one for each triangle sharing it.
994
995            If there are multiple cells with same maximum value, the
996            first cell encountered in the triangle array is returned.
997        """
998
999        i = self.get_maximum_index(indices)
1000        x, y = self.domain.get_centroid_coordinates()[i]
1001
1002        return x, y
1003
1004    ##
1005    # @brief  Get index for minimum value of quantity.
1006    # @param indices Set of IDs of elements to work on.
1007    def get_minimum_index(self, indices=None):
1008        """See get extreme index for details"""
1009
1010        return self.get_extremum_index(mode='min', indices=indices)
1011
1012    ##
1013    # @brief Return minimum value of quantity (on centroids).
1014    # @param indices Set of IDs of elements to work on.
1015    def get_minimum_value(self, indices=None):
1016        """Return minimum value of quantity (on centroids)
1017
1018        Optional argument:
1019            indices is the set of element ids that the operation applies to.
1020
1021        Usage:
1022            v = get_minimum_value()
1023
1024        See get_maximum_value for more details.
1025        """
1026
1027        i = self.get_minimum_index(indices)
1028        V = self.get_values(location='centroids')
1029
1030        return V[i]
1031
1032
1033    ##
1034    # @brief Get location of minimum value of quantity (on centroids).
1035    # @param indices Set of IDs of elements to work on.
1036    def get_minimum_location(self, indices=None):
1037        """Return location of minimum value of quantity (on centroids)
1038
1039        Optional argument:
1040            indices is the set of element ids that the operation applies to.
1041
1042        Usage:
1043            x, y = get_minimum_location()
1044
1045        Notes:
1046            We do not seek the maximum at vertices as each vertex can
1047            have multiple values - one for each triangle sharing it.
1048
1049            If there are multiple cells with same maximum value, the
1050            first cell encountered in the triangle array is returned.
1051        """
1052
1053        i = self.get_minimum_index(indices)
1054        x, y = self.domain.get_centroid_coordinates()[i]
1055
1056        return x, y
1057
1058    ##
1059    # @brief Get values at interpolation points.
1060    # @param interpolation_points List of UTM coords or geospatial data object.
1061    # @param use_cache ??
1062    # @param verbose True if this method is to be verbose.
1063    def get_interpolated_values(self, interpolation_points,
1064                                      use_cache=False,
1065                                      verbose=False):
1066        """Get values at interpolation points
1067
1068        The argument interpolation points must be given as either a
1069        list of absolute UTM coordinates or a geospatial data object.
1070        """
1071
1072        # FIXME (Ole): Points might be converted to coordinates relative to mesh origin
1073        # This could all be refactored using the
1074        # 'change_points_geo_ref' method of Class geo_reference.
1075        # The purpose is to make interpolation points relative
1076        # to the mesh origin.
1077        #
1078        # Speed is also a consideration here.
1079
1080        # Ensure that interpolation points is either a list of
1081        # points, Nx2 array, or geospatial and convert to Numeric array
1082        if isinstance(interpolation_points, Geospatial_data):
1083            # Ensure interpolation points are in absolute UTM coordinates
1084            interpolation_points = \
1085                interpolation_points.get_data_points(absolute=True)
1086
1087        # Reconcile interpolation points with georeference of domain
1088        interpolation_points = \
1089            self.domain.geo_reference.get_relative(interpolation_points)
1090        interpolation_points = ensure_numeric(interpolation_points)
1091
1092
1093        # Get internal representation (disconnected) of vertex values
1094        vertex_values, triangles = self.get_vertex_values(xy=False,
1095                                                          smooth=False)
1096
1097        # Get possibly precomputed interpolation object
1098        I = self.domain.get_interpolation_object()
1099
1100        # Call interpolate method with interpolation points
1101        result = I.interpolate_block(vertex_values, interpolation_points,
1102                                     use_cache=use_cache, verbose=verbose)
1103
1104        return result
1105
1106    ##
1107    # @brief Get values as an array.
1108    # @param interpolation_points List of coords to get values at.
1109    # @param location Where to store results.
1110    # @param indices Set of IDs of elements to work on.
1111    # @param use_cache
1112    # @param verbose True if this method is to be verbose.
1113    def get_values(self, 
1114                   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        # YES (Ole): Edge values are necessary for volumetric balance check and inflow boundary. Keep them.
1172
1173        if location not in ['vertices', 'centroids',
1174                            'edges', 'unique vertices']:
1175            msg = 'Invalid location: %s' % location
1176            raise msg
1177
1178        import types
1179
1180        assert type(indices) in [types.ListType, types.NoneType,
1181                                 num.ArrayType], \
1182                   'Indices must be a list or None'
1183
1184        if location == 'centroids':
1185            if (indices ==  None):
1186                indices = range(len(self))
1187            return num.take(self.centroid_values, indices)
1188        elif location == 'edges':
1189            if (indices ==  None):
1190                indices = range(len(self))
1191            return num.take(self.edge_values, indices)
1192        elif location == 'unique vertices':
1193            if (indices ==  None):
1194                indices=range(self.domain.get_number_of_nodes())
1195            vert_values = []
1196
1197            # Go through list of unique vertices
1198            for unique_vert_id in indices:
1199                triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1200
1201                # In case there are unused points
1202                if len(triangles) == 0:
1203                    msg = 'Unique vertex not associated with triangles'
1204                    raise Exception, msg
1205
1206                # Go through all triangle, vertex pairs
1207                # Average the values
1208                # FIXME (Ole): Should we merge this with get_vertex_values
1209                sum = 0
1210                for triangle_id, vertex_id in triangles:
1211                    sum += self.vertex_values[triangle_id, vertex_id]
1212                vert_values.append(sum / len(triangles))
1213            return num.array(vert_values, num.Float)
1214        else:
1215            if (indices is None):
1216                indices = range(len(self))
1217            return num.take(self.vertex_values, indices)
1218
1219    ##
1220    # @brief Set vertex values for all unique vertices based on array.
1221    # @param A Array to set values with.
1222    # @param indices Set of IDs of elements to work on.
1223    # @param use_cache ??
1224    # @param verbose??
1225    def set_vertex_values(self, A,
1226                                indices=None,
1227                                use_cache=False,
1228                                verbose=False):
1229        """Set vertex values for all unique vertices based on input array A
1230        which has one entry per unique vertex, i.e. one value for each row in
1231        array self.domain.nodes.
1232
1233        indices is the list of vertex_id's that will be set.
1234
1235        This function is used by set_values_from_array
1236        """
1237
1238        # Check that A can be converted to array and is of appropriate dim
1239        A = ensure_numeric(A, num.Float)
1240        assert len(A.shape) == 1
1241
1242        if indices is None:
1243            assert A.shape[0] == self.domain.get_nodes().shape[0]
1244            vertex_list = range(A.shape[0])
1245        else:
1246            assert A.shape[0] == len(indices)
1247            vertex_list = indices
1248
1249        #FIXME(Ole): This function ought to be faster.
1250        # We need to get the triangles_and_vertices list
1251        # from domain in one hit, then cache the computation of the
1252        # Nx3 array of vertex values that can then be assigned using
1253        # set_values_from_array.
1254        #
1255        # Alternatively, some C code would be handy
1256        #
1257        self._set_vertex_values(vertex_list, A)
1258           
1259    ##
1260    # @brief Go through list of unique vertices.
1261    # @param vertex_list ??
1262    # @param A ??
1263    def _set_vertex_values(self, vertex_list, A):
1264        """Go through list of unique vertices
1265        This is the common case e.g. when values
1266        are obtained from a pts file through fitting
1267        """
1268
1269        # Go through list of unique vertices
1270        for i_index, unique_vert_id in enumerate(vertex_list):
1271            triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1272
1273            # In case there are unused points
1274            if len(triangles) == 0:
1275                continue
1276
1277            # Go through all triangle, vertex pairs
1278            # touching vertex unique_vert_id and set corresponding vertex value
1279            for triangle_id, vertex_id in triangles:
1280                self.vertex_values[triangle_id, vertex_id] = A[i_index]
1281
1282        # Intialise centroid and edge_values
1283        self.interpolate()
1284
1285    ##
1286    # @brief Smooth vertex values.
1287    def smooth_vertex_values(self, use_cache=False, verbose=False):
1288        """Smooths vertex values."""
1289
1290        A, V = self.get_vertex_values(xy=False, smooth=True)
1291        self.set_vertex_values(A, use_cache=use_cache, verbose=verbose)
1292
1293    ############################################################################
1294    # Methods for outputting model results
1295    ############################################################################
1296
1297    ##
1298    # @brief Get vertex values like an OBJ format i.e. one value per node.
1299    # @param xy True if we return X and Y as well as A and V.
1300    # @param smooth True if vertex values are to be smoothed.
1301    # @param precision The type of the result values (default float).
1302    def get_vertex_values(self, xy=True, smooth=None, precision=None):
1303        """Return vertex values like an OBJ format i.e. one value per node.
1304
1305        The vertex values are returned as one sequence in the 1D float array A.
1306        If requested the coordinates will be returned in 1D arrays X and Y.
1307
1308        The connectivity is represented as an integer array, V, of dimension
1309        Mx3, where M is the number of triangles. Each row has three indices
1310        defining the triangle and they correspond to elements in the arrays
1311        X, Y and A.
1312
1313        If smooth is True, vertex values corresponding to one common coordinate
1314        set will be smoothed by taking the average of vertex values for each
1315        node.  In this case vertex coordinates will be de-duplicated
1316        corresponding to the original nodes as obtained from the method
1317        general_mesh.get_nodes()
1318
1319        If no smoothings is required, vertex coordinates and values will be
1320        aggregated as a concatenation of values at vertices 0, vertices 1 and
1321        vertices 2.  This corresponds to the node coordinates obtained from the
1322        method general_mesh.get_vertex_coordinates()
1323
1324        Calling convention
1325        if xy is True:
1326           X, Y, A, V = get_vertex_values
1327        else:
1328           A, V = get_vertex_values
1329        """
1330
1331        if smooth is None:
1332            # Take default from domain
1333            try:
1334                smooth = self.domain.smooth
1335            except:
1336                smooth = False
1337
1338        if precision is None:
1339            precision = num.Float
1340
1341        if smooth is True:
1342            # Ensure continuous vertex values by averaging values at each node
1343            V = self.domain.get_triangles()
1344            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1345            A = num.zeros(N, num.Float)
1346            points = self.domain.get_nodes()
1347
1348            if 1:
1349                # Fast C version
1350                average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1351                                      ensure_numeric(self.domain.number_of_triangles_per_node),
1352                                      ensure_numeric(self.vertex_values),
1353                                      A)
1354                A = A.astype(precision)
1355            else:
1356                # Slow Python version
1357                current_node = 0
1358                k = 0 # Track triangles touching on node
1359                total = 0.0
1360                for index in self.domain.vertex_value_indices:
1361                    if current_node == N:
1362                        msg = 'Current node exceeding number of nodes (%d) ' % N
1363                        raise msg
1364
1365                    k += 1
1366
1367                    volume_id = index / 3
1368                    vertex_id = index % 3
1369
1370                    v = self.vertex_values[volume_id, vertex_id]
1371                    total += v
1372
1373                    if self.domain.number_of_triangles_per_node[current_node] == k:
1374                        A[current_node] = total/k
1375
1376                        # Move on to next node
1377                        total = 0.0
1378                        k = 0
1379                        current_node += 1
1380        else:
1381            # Return disconnected internal vertex values
1382            V = self.domain.get_disconnected_triangles()
1383            points = self.domain.get_vertex_coordinates()
1384            A = self.vertex_values.flat.astype(precision)
1385
1386        # Return
1387        if xy is True:
1388            X = points[:,0].astype(precision)
1389            Y = points[:,1].astype(precision)
1390
1391            return X, Y, A, V
1392        else:
1393            return A, V
1394
1395    ##
1396    # @brief Extrapolate conserved quantities from centroid.
1397    def extrapolate_first_order(self):
1398        """Extrapolate conserved quantities from centroid to vertices and edges
1399        for each volume using first order scheme.
1400        """
1401
1402        qc = self.centroid_values
1403        qv = self.vertex_values
1404        qe = self.edge_values
1405
1406        for i in range(3):
1407            qv[:,i] = qc
1408            qe[:,i] = qc
1409
1410        self.x_gradient *= 0.0
1411        self.y_gradient *= 0.0
1412
1413    ##
1414    # @brief Compute the integral of quantity across entire domain.
1415    # @return The integral.
1416    def get_integral(self):
1417        """Compute the integral of quantity across entire domain."""
1418
1419        areas = self.domain.get_areas()
1420        integral = 0
1421        for k in range(len(self.domain)):
1422            area = areas[k]
1423            qc = self.centroid_values[k]
1424            integral += qc*area
1425
1426        return integral
1427
1428    ##
1429    # @brief get the gradients.
1430    def get_gradients(self):
1431        """Provide gradients. Use compute_gradients first."""
1432
1433        return self.x_gradient, self.y_gradient
1434
1435    ##
1436    # @brief ??
1437    # @param timestep ??
1438    def update(self, timestep):
1439        # Call correct module function
1440        # (either from this module or C-extension)
1441        return update(self, timestep)
1442
1443    ##
1444    # @brief ??
1445    def compute_gradients(self):
1446        # Call correct module function
1447        # (either from this module or C-extension)
1448        return compute_gradients(self)
1449
1450    ##
1451    # @brief ??
1452    def limit(self):
1453        # Call correct module depending on whether
1454        # basing limit calculations on edges or vertices
1455        limit_old(self)
1456
1457    ##
1458    # @brief ??
1459    def limit_vertices_by_all_neighbours(self):
1460        # Call correct module function
1461        # (either from this module or C-extension)
1462        limit_vertices_by_all_neighbours(self)
1463
1464    ##
1465    # @brief ??
1466    def limit_edges_by_all_neighbours(self):
1467        # Call correct module function
1468        # (either from this module or C-extension)
1469        limit_edges_by_all_neighbours(self)
1470
1471    ##
1472    # @brief ??
1473    def limit_edges_by_neighbour(self):
1474        # Call correct module function
1475        # (either from this module or C-extension)
1476        limit_edges_by_neighbour(self)
1477
1478    ##
1479    # @brief ??
1480    def extrapolate_second_order(self):
1481        # Call correct module function
1482        # (either from this module or C-extension)
1483        compute_gradients(self)
1484        extrapolate_from_gradient(self)
1485
1486    ##
1487    # @brief ??
1488    def extrapolate_second_order_and_limit_by_edge(self):
1489        # Call correct module function
1490        # (either from this module or C-extension)
1491        extrapolate_second_order_and_limit_by_edge(self)
1492
1493    ##
1494    # @brief ??
1495    def extrapolate_second_order_and_limit_by_vertex(self):
1496        # Call correct module function
1497        # (either from this module or C-extension)
1498        extrapolate_second_order_and_limit_by_vertex(self)
1499
1500    ##
1501    # @brief ??
1502    # @param bound ??
1503    def bound_vertices_below_by_constant(self, bound):
1504        # Call correct module function
1505        # (either from this module or C-extension)
1506        bound_vertices_below_by_constant(self, bound)
1507
1508    ##
1509    # @brief ??
1510    # @param quantity ??
1511    def bound_vertices_below_by_quantity(self, quantity):
1512        # Call correct module function
1513        # (either from this module or C-extension)
1514
1515        # check consistency
1516        assert self.domain == quantity.domain
1517        bound_vertices_below_by_quantity(self, quantity)
1518
1519    ##
1520    # @brief ??
1521    def backup_centroid_values(self):
1522        # Call correct module function
1523        # (either from this module or C-extension)
1524        backup_centroid_values(self)
1525
1526    ##
1527    # @brief ??
1528    # @param a ??
1529    # @param b ??
1530    def saxpy_centroid_values(self, a, b):
1531        # Call correct module function
1532        # (either from this module or C-extension)
1533        saxpy_centroid_values(self, a, b)
1534
1535
1536##
1537# @brief OBSOLETE!
1538class Conserved_quantity(Quantity):
1539    """Class conserved quantity being removed, use Quantity."""
1540
1541    def __init__(self, domain, vertex_values=None):
1542        msg = 'ERROR: Use Quantity instead of Conserved_quantity'
1543        raise Exception, msg
1544
1545
1546######
1547# Prepare the C extensions.
1548######
1549
1550from anuga.utilities import compile
1551
1552if compile.can_use_C_extension('quantity_ext.c'):
1553    # Underlying C implementations can be accessed
1554
1555    from quantity_ext import \
1556         average_vertex_values,\
1557         backup_centroid_values,\
1558         saxpy_centroid_values,\
1559         compute_gradients,\
1560         limit_old,\
1561         limit_vertices_by_all_neighbours,\
1562         limit_edges_by_all_neighbours,\
1563         limit_edges_by_neighbour,\
1564         limit_gradient_by_neighbour,\
1565         extrapolate_from_gradient,\
1566         extrapolate_second_order_and_limit_by_edge,\
1567         extrapolate_second_order_and_limit_by_vertex,\
1568         bound_vertices_below_by_constant,\
1569         bound_vertices_below_by_quantity,\
1570         interpolate_from_vertices_to_edges,\
1571         interpolate_from_edges_to_vertices,\
1572         update
1573else:
1574    msg = 'C implementations could not be accessed by %s.\n ' % __file__
1575    msg += 'Make sure compile_all.py has been run as described in '
1576    msg += 'the ANUGA installation guide.'
1577    raise Exception, msg
Note: See TracBrowser for help on using the repository browser.