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

Last change on this file since 4886 was 4886, checked in by steve, 16 years ago

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