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

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

After changes to get_absolute, ensure_numeric, etc.

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