Changeset 279 for inundation/ga


Ignore:
Timestamp:
Sep 7, 2004, 5:24:43 PM (20 years ago)
Author:
steve
Message:

Working on Quantities

Location:
inundation/ga/storm_surge/pyvolution-1d
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution-1d/domain.py

    r256 r279  
    2626
    2727        #Allocate space for geometric quantities
    28         self.vertex_coordinates   = zeros((N, 2), Float)
    29         self.centroid_coordinates = zeros(N, Float)
    30         self.areas                = zeros(N, Float)
     28        self.vertices  = zeros((N, 2), Float)
     29        self.centroids = zeros(N, Float)
     30        self.areas     = zeros(N, Float)
    3131        for i in range(N):
    3232            xl = self.coordinates[i]
    3333            xr = self.coordinates[i+1]
    34             self.vertex_coordinates[i,0] = xl
    35             self.vertex_coordinates[i,1] = xr
     34            self.vertices[i,0] = xl
     35            self.vertices[i,1] = xr
    3636           
    3737            centroid = (xl+xr)/2
    38             self.centroid_coordinates[i] = centroid
     38            self.centroids[i] = centroid
    3939
    4040            msg = 'Coordinates should be ordered, smallest to largest'
    … …  
    4343            self.areas[i] = (xr-xl)
    4444
     45##         print 'N', N
     46##         print 'Centroid', self.centroids
     47##         print 'Areas', self.areas
     48##         print 'Vertex_Coordinates', self.vertices
     49       
     50    def get_centroids(self):
     51        """Return all coordinates of centroids
     52        Return x coordinate of centroid for each element as a N array
     53        """
     54       
     55        return self.centroids
    4556
    46         print 'N', N
    47         print 'Coordinates', self.coordinates
    48         print 'Centroid', self.centroid_coordinates
    49         print 'Areas', self.areas
    50         print 'Vertex_Coordinates', self.vertex_coordinates
    51        
    52     def get_centroid_coordinates(self):
     57    def get_vertices(self):
    5358        """Return all coordinates of centroids
    5459        Return x coordinate of centroid for each element as a N array
    5560        """
    5661
    57         return self.centroid_coordinates
     62        return self.vertices
    5863
    5964    def get_coordinate(self, elem_id, vertex=None):
    … …  
    6469
    6570        if vertex is None:
    66             return self.centroid_coordinates[elem_id]
     71            return self.centroids[elem_id]
    6772        else:
    68             return self.vertex_coordinates[elem_id,vertex]
    69        
     73            return self.vertices[elem_id,vertex]
     74
     75    def get_area(self, elem_id):
     76        """Return area of element id
     77        """
     78
     79        return self.areas[elem_id]
     80   
    7081
    7182if __name__ == "__main__":

  • inundation/ga/storm_surge/pyvolution-1d/quantity.py

    r256 r279  
    5050        self.centroid_values = zeros(N, Float)
    5151
    52         #Intialise centroid and edge_values
     52        #Intialise centroid values
    5353        self.interpolate()
    54 
    55 
    56 
     54       
    5755   
    5856    def interpolate(self):
    … …  
    6866            self.centroid_values[i] = (v0 + v1)/2
    6967
    70 if __name__ == "__main__":
    71     from domain import Domain
    72 
    73     points1 = [0.0, 1.0, 2.0, 3.0]
    74     vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
    75    
    76     D1 = Domain(points1)
    77 
    78     Q1 = Quantity(D1, vertex_values)
    79 
    80     print Q1.vertex_values
    81     print Q1.centroid_values
    82    
    83 
    84 
    85            
    86            
    8768    def set_values(self, X, location='vertices'):
    8869        """Set values for quantity
    … …  
    9576        In case of location == 'centroid' the dimension values must
    9677        be a list of a Numerical array of length N, N being the number
    97         of elements in the mesh. Otherwise it must be of dimension Nx2
     78        of elements in the mesh. Otherwise it must be of dimension Nx3
    9879
    9980        The values will be stored in elements following their
    … …  
    10283        If values are described a function, it will be evaluated at specified points
    10384
    104         If selected location is vertices, values for centroid
     85        If selected location is vertices, values for centroid and edges
    10586        will be assigned interpolated values.
    10687        In any other case, only values for the specified locations
    … …  
    10990
    11091        if location not in ['vertices', 'centroids']:
    111             msg = 'Invalid location: %s' %location
     92            msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
    11293            raise msg
    11394
    … …  
    126107            else:
    127108                self.vertex_values[:] = X               
     109
    128110        else:
    129111            #Use array specific method
    … …  
    134116            self.interpolate()
    135117           
     118
    136119
    137120
    … …  
    150133        else:
    151134            #Vertices
    152             P = self.domain.get_vertex_coordinates()
     135            P = self.domain.get_vertices()
    153136            for i in range(2):
    154                 self.vertex_values[:,i] = f(P[:,2*i])
    155 
    156            
    157 ##     def set_array_values(self, values, location='vertices'):
    158 ##         """Set values for quantity
    159 
    160 ##         values: Numeric array
    161 ##         location: Where values are to be stored.
    162 ##                   Permissible options are: vertices, edges, centroid
    163 ##                   Default is "vertices"
    164 
    165 ##         In case of location == 'centroid' the dimension values must
    166 ##         be a list of a Numerical array of length N, N being the number
    167 ##         of elements in the mesh. Otherwise it must be of dimension Nx3
    168 
    169 ##         The values will be stored in elements following their
    170 ##         internal ordering.
    171 
    172 ##         If selected location is vertices, values for centroid and edges
    173 ##         will be assigned interpolated values.
    174 ##         In any other case, only values for the specified locations
    175 ##         will be assigned and the others will be left undefined.
    176 ##         """
    177 
    178 ##         from Numeric import array, Float
    179 
    180 ##         values = array(values).astype(Float)
    181 
    182 ##         N = self.centroid_values.shape[0]
    183        
    184 ##         msg = 'Number of values must match number of elements'
    185 ##         assert values.shape[0] == N, msg
    186        
    187 ##         if location == 'centroids':
    188 ##             assert len(values.shape) == 1, 'Values array must be 1d'
    189 ##             self.centroid_values = values                       
    190 ##         elif location == 'edges':
    191 ##             assert len(values.shape) == 2, 'Values array must be 2d'
    192 ##             msg = 'Array must be N x 3'           
    193 ##             assert values.shape[1] == 3, msg
    194            
    195 ##             self.edge_values = values
    196 ##         else:
    197 ##             assert len(values.shape) == 2, 'Values array must be 2d'
    198 ##             msg = 'Array must be N x 3'           
    199 ##             assert values.shape[1] == 3, msg
    200            
    201 ##             self.vertex_values = values
    202 
     137                self.vertex_values[:,i] = f(P[:,i])
     138
     139
     140    def set_array_values(self, values, location='vertices'):
     141        """Set values for quantity
     142
     143        values: Numeric array
     144        location: Where values are to be stored.
     145                  Permissible options are: vertices, centroid
     146                  Default is "vertices"
     147
     148        In case of location == 'centroid' the dimension values must
     149        be a list of a Numerical array of length N, N being the number
     150        of elements in the mesh. Otherwise it must be of dimension Nx2
     151
     152        The values will be stored in elements following their
     153        internal ordering.
     154
     155        If selected location is vertices, values for centroid
     156        will be assigned interpolated values.
     157        In any other case, only values for the specified locations
     158        will be assigned and the others will be left undefined.
     159        """
     160
     161        from Numeric import array, Float
     162
     163        values = array(values).astype(Float)
     164
     165        N = self.centroid_values.shape[0]
     166       
     167        msg = 'Number of values must match number of elements'
     168        assert values.shape[0] == N, msg
     169       
     170        if location == 'centroids':
     171            assert len(values.shape) == 1, 'Values array must be 1d'
     172            self.centroid_values = values                       
     173        else:
     174            assert len(values.shape) == 2, 'Values array must be 2d'
     175            msg = 'Array must be N x 2'           
     176            assert values.shape[1] == 2, msg
     177           
     178            self.vertex_values = values
     179
     180
     181
     182if __name__ == "__main__":
     183    from domain import Domain
     184
     185    points1 = [0.0, 1.0, 2.0, 3.0]
     186    vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
     187   
     188    D1 = Domain(points1)
     189
     190    Q1 = Quantity(D1, vertex_values)
     191
     192    print Q1.vertex_values
     193    print Q1.centroid_values
     194
     195    new_vertex_values = [[2.0,1.0],[3.0,4.0],[-2.0,4.0]]
     196
     197    Q1.set_values(new_vertex_values)
     198   
     199   
     200    print Q1.vertex_values
     201    print Q1.centroid_values
     202
     203    new_centroid_values = [20,30,40]
     204    Q1.set_values(new_centroid_values,'centroids')
     205           
     206    print Q1.vertex_values
     207    print Q1.centroid_values
     208
     209    def fun(x):
     210
     211        return x**2
     212
     213    Q1.set_values(fun,'vertices')   
     214           
     215    print Q1.vertex_values
     216    print Q1.centroid_values
    203217
    204218

  • inundation/ga/storm_surge/pyvolution-1d/test_domain.py

    r256 r279  
    1010class TestCase(unittest.TestCase):
    1111    def setUp(self):
    12        
     12        pass
     13
     14
     15    def tearDown(self):
     16        pass
     17
     18
     19    def test_construct(self):
     20
    1321        a = 0.0
    1422        b = 2.0
    … …  
    2028        points = [a, b, c, d, e, f]
    2129
    22         self.points = [a, b, c, d, e, f]
    23         self.domain = Domain(self.points)
    24 
    25        
    26     def tearDown(self):
    27         pass
     30        domain = Domain(points)
     31        assert allclose(domain.get_centroids(),[(a+b)*0.5, (b+c)*0.5, (c+d)*0.5, (d+e)*0.5, (e+f)*0.5])
    2832
    2933
    30     def test_construct(self):
     34    def test_area(self):
    3135
    32         domain = self.domain
    33         centroids = domain.get_centroids()
    34         assert allclose(centroids, [1.0, 2.25, 2.75, 6.5, 11.0])
    35            
     36        a = 0.0
     37        b = 2.0
     38        c = 2.5
     39        d = 3.0
     40        e = 10.0
     41        f = 12.0
     42
     43        points = [a, b, c, d, e, f]
     44
     45        domain = Domain(points)
     46        assert allclose(domain.get_area(0),b-a)
     47        assert allclose(domain.get_area(1),c-b)
     48        assert allclose(domain.get_area(2),d-c)
     49        assert allclose(domain.get_area(3),e-d)
     50        assert allclose(domain.get_area(4),f-e)
     51       
    3652    def test_in_order(self):
    3753        a = 0.0
    … …  
    6278        domain = Domain(points)
    6379
    64         assert allclose(domain.get_point(0), 0.5*(a+b))
    65         assert allclose(domain.get_point(1), 0.5*(b+c))
    66         assert allclose(domain.get_point(2), 0.5*(c+d))
     80        assert allclose(domain.get_coordinate(0), 0.5*(a+b))
     81        assert allclose(domain.get_coordinate(1), 0.5*(b+c))
     82        assert allclose(domain.get_coordinate(2), 0.5*(c+d))
    6783       
    6884        try:
    69             allclose(domain.get_point(3), 0.5*(a+b))
     85            allclose(domain.get_coordinate(3), 0.5*(a+b))
    7086        except:
    7187            pass
    … …  
    7490            raise msg
    7591       
    76         assert allclose(domain.get_point(0,0), a)
    77         assert allclose(domain.get_point(1,0), b)
    78         assert allclose(domain.get_point(2,0), c)
     92        assert allclose(domain.get_coordinate(0,0), a)
     93        assert allclose(domain.get_coordinate(1,0), b)
     94        assert allclose(domain.get_coordinate(2,0), c)
    7995
    80         assert allclose(domain.get_point(0,1), b)
    81         assert allclose(domain.get_point(1,1), c)
    82         assert allclose(domain.get_point(2,1), d)
     96        assert allclose(domain.get_coordinate(0,1), b)
     97        assert allclose(domain.get_coordinate(1,1), c)
     98        assert allclose(domain.get_coordinate(2,1), d)
    8399
    84100

  • inundation/ga/storm_surge/pyvolution-1d/test_quantity.py

    r256 r279  
    8585
    8686
    87 ##     def test_set_values(self):
    88 ##         quantity = Quantity(self.mesh4)
    89 
    90 
    91 ##         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
    92 ##                             location = 'vertices')
    93 ##         assert allclose(quantity.vertex_values,
    94 ##                         [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    95 ##         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    96 ##         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    97 ##                                                [5., 5., 5.],
    98 ##                                                [4.5, 4.5, 0.],
    99 ##                                                [3.0, -1.5, -1.5]])
    100 
    101 
    102 ##         #Test default
    103 ##         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    104 ##         assert allclose(quantity.vertex_values,
    105 ##                         [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
    106 ##         assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    107 ##         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    108 ##                                                [5., 5., 5.],
    109 ##                                                [4.5, 4.5, 0.],
    110 ##                                                [3.0, -1.5, -1.5]])
    111 
    112 ##         #Test centroids
    113 ##         quantity.set_values([1,2,3,4], location = 'centroids')
    114 ##         assert allclose(quantity.centroid_values, [1., 2., 3., 4.]) #Centroid
    115 
    116 ##         #Test edges
    117 ##         quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
    118 ##                             location = 'edges')
    119 ##         assert allclose(quantity.edge_values,
    120 ##                         [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])       
    121 
    122 ##         #Test exceptions
    123 ##         try:
    124 ##             quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
    125 ##                                 location = 'bas kamel tuba')
    126 ##         except:
    127 ##             pass
    128 
    129 
    130 ##         try:
    131 ##             quantity.set_values([[1,2,3], [0,0,9]])
    132 ##         except AssertionError:
    133 ##             pass
    134 ##         except:
    135 ##             raise 'should have raised Assertionerror'
    136 
    137 
    138 
    139 ##     def test_set_values_const(self):
    140 ##         quantity = Quantity(self.mesh4)
    141 
    142 ##         quantity.set_values(1.0, location = 'vertices')
    143 ##         assert allclose(quantity.vertex_values,
    144 ##                         [[1,1,1], [1,1,1], [1,1,1], [1, 1, 1]])
    145 ##         assert allclose(quantity.centroid_values, [1, 1, 1, 1]) #Centroid
    146 ##         assert allclose(quantity.edge_values, [[1, 1, 1],
    147 ##                                                [1, 1, 1],
    148 ##                                                [1, 1, 1],
    149 ##                                                [1, 1, 1]])
    150 
    151 
    152 ##         quantity.set_values(2.0, location = 'centroids')
    153 ##         assert allclose(quantity.centroid_values, [2, 2, 2, 2])
    154 
    155 ##         quantity.set_values(3.0, location = 'edges')
    156 ##         assert allclose(quantity.edge_values, [[3, 3, 3],
    157 ##                                                [3, 3, 3],
    158 ##                                                [3, 3, 3],
    159 ##                                                [3, 3, 3]])       
    160 
    161 
    162 ##     def test_set_values_func(self):
    163 ##         quantity = Quantity(self.mesh4)
    164 
    165 ##         def f(x, y):
    166 ##             return x+y
    167 
    168 ##         quantity.set_values(f, location = 'vertices')
    169 ##         assert allclose(quantity.vertex_values,
    170 ##                         [[2,0,2], [2,2,4], [4,2,4], [4,2,4]])       
    171 ##         assert allclose(quantity.centroid_values,
    172 ##                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])
    173 ##         assert allclose(quantity.edge_values,
    174 ##                         [[1,2,1], [3,3,2], [3,4,3], [3,4,3]])               
    175 
    176        
    177 ##         quantity.set_values(f, location = 'centroids')
    178 ##         assert allclose(quantity.centroid_values,
    179 ##                         [4.0/3, 8.0/3, 10.0/3, 10.0/3])       
     87    def test_set_values(self):
     88        quantity = Quantity(self.domain5)
     89
     90
     91        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
     92                            location = 'vertices')
     93        assert allclose(quantity.vertex_values,
     94                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
     95        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
     96
     97        #Test default
     98        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
     99        assert allclose(quantity.vertex_values,
     100                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])       
     101        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid
     102
     103        #Test centroids
     104        quantity.set_values([1,2,3,4,5], location = 'centroids')
     105        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid
     106
     107        #Test exceptions
     108        try:
     109            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
     110                                location = 'bas kamel tuba')
     111        except:
     112            pass
     113
     114
     115        try:
     116            quantity.set_values([[1,2], [0,0]])
     117        except AssertionError:
     118            pass
     119        except:
     120            raise 'should have raised AssertionError'
     121
     122
     123
     124    def test_set_values_const(self):
     125        quantity = Quantity(self.domain5)
     126
     127        quantity.set_values(1.0, location = 'vertices')
     128        assert allclose(quantity.vertex_values,
     129                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
     130        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid
     131
     132
     133        quantity.set_values(2.0, location = 'centroids')
     134        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])
     135
     136
     137    def test_set_values_func(self):
     138        quantity = Quantity(self.domain5)
     139
     140        def f(x):
     141            return x**2
     142
     143        quantity.set_values(f, location = 'vertices')
     144        assert allclose(quantity.vertex_values,
     145                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])       
     146        assert allclose(quantity.centroid_values,
     147                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
     148       
     149        quantity.set_values(f, location = 'centroids')
     150        assert allclose(quantity.centroid_values,
     151                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])       
    180152
    181153
    … …  
    223195##         assert b[0] == 0.0
    224196##         #The others are OK
    225 ##         for i in range(1,4):
     197##      163.968025   for i in range(1,4):
    226198##             assert allclose(a[i], 3.0)
    227199##             assert allclose(b[i], 1.0)
    … …  
    303275## #         assert v2.conserved_quantities_vertex2 <= qmax
    304276## #         assert v2.conserved_quantities_vertex2 >= qmin                   
    305 
     277163.968025
    306278
    307279## #         #Check that volume has been preserved   
Note: See TracChangeset for help on using the changeset viewer.