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

Last change on this file since 4839 was 4839, checked in by duncan, 16 years ago

speeding up slow functions in general mesh. Fix for ticket 201. Also checking in GA validation test for profiling fitting

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