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

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

Refactored Mesh-Domain inheritance pattern into a composition pattern, thereby
paving the way for reuse of Mesh instance in fitting as per ticket:242.
Had to disable test for isinstance(domain, Domain) in quantity.py as it
failed for unknown reasons. All other tests and validation suite passes.

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