# source:inundation/pyvolution/quantity.py@1903

Last change on this file since 1903 was 1903, checked in by duncan, 18 years ago

removing / gutting duplication of methods reading xya files

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