source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py @ 8123

Last change on this file since 8123 was 8123, checked in by steve, 13 years ago

Added in a procedure to calculate local gradients on each triangle using
vectex values

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