source: inundation/pyvolution/quantity.py @ 2648

Last change on this file since 2648 was 2648, checked in by steve, 18 years ago

Looking at island.py example. I changed the order of the implicit and explicit update (in quantity_ext.c) which I think improved the situation a little. But this has left a few fails in the test_all suite which we need to fix up. Mainly small differences in results.

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