Changeset 8930


Ignore:
Timestamp:
Jun 26, 2013, 9:29:55 AM (12 years ago)
Author:
steve
Message:

Updated versoin of set_elevation etc

Location:
trunk/anuga_core/source
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r8871 r8930  
    239239            Q.set_values(other)
    240240
    241         result = Quantity(self.domain)
    242241
    243242        # The maximum of vertex_values, edge_values and centroid_values
     
    245244        # set_values (which calls interpolate). Otherwise
    246245        # edge and centroid values wouldn't be max from q1 and q2
    247         result.vertex_values = num.maximum(self.vertex_values, Q.vertex_values)
    248         result.edge_values = num.maximum(self.edge_values, Q.edge_values)
    249         result.centroid_values = num.maximum(self.centroid_values, Q.centroid_values)
    250 
    251         return result
     246        self.vertex_values[:] = num.maximum(self.vertex_values, Q.vertex_values)
     247        self.edge_values[:] = num.maximum(self.edge_values, Q.edge_values)
     248        self.centroid_values[:] = num.maximum(self.centroid_values, Q.centroid_values)
     249
     250       
    252251
    253252
     
    266265            Q.set_values(other)
    267266
    268         result = Quantity(self.domain)
    269267
    270268        # The minimum of vertex_values, edge_values and centroid_values
    271269        # are calculated and assigned directly without using
    272270        # set_values (which calls interpolate). Otherwise
    273         result.vertex_values = num.minimum(self.vertex_values, Q.vertex_values)
    274         result.edge_values = num.minimum(self.edge_values, Q.edge_values)
    275         result.centroid_values = num.minimum(self.centroid_values, Q.centroid_values)
    276 
    277         return result
     271        self.vertex_values[:]= num.minimum(self.vertex_values, Q.vertex_values)
     272        self.edge_values[:] = num.minimum(self.edge_values, Q.edge_values)
     273        self.centroid_values[:] = num.minimum(self.centroid_values, Q.centroid_values)
     274
     275
    278276
    279277
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r8486 r8930  
    24382438                           
    24392439
     2440    def test_maximum(self):
     2441        quantity = Quantity(self.mesh4)
     2442
     2443        # get referece to data arrays
     2444        centroid_values = quantity.centroid_values
     2445        vertex_values = quantity.vertex_values
     2446        edge_values = quantity.edge_values
     2447
     2448        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
     2449                            location = 'vertices')
     2450        assert num.allclose(quantity.vertex_values,
     2451                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
     2452
     2453        assert id(vertex_values) == id(quantity.vertex_values)
     2454        assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
     2455        assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
     2456                                                   [5., 5., 5.],
     2457                                                   [4.5, 4.5, 0.],
     2458                                                   [3.0, -1.5, -1.5]])
     2459
     2460
     2461
     2462        other_quantity = Quantity(self.mesh4)
     2463
     2464        other_quantity.set_values([[0,0,0], [1,1,6], [10,10,10], [0,0,4]],
     2465                            location = 'vertices')
     2466
     2467 
     2468        #===============================
     2469        quantity.maximum(other_quantity)
     2470        #===============================
     2471
     2472        exact_vertex_values = num.array([[  1.,   2.,   3.],
     2473                                         [  5.,   5.,  6.],
     2474                                         [ 10.,  10.,  10.],
     2475                                         [  0.,   3.,  4.]])
     2476        exact_centroid_values = num.array([  2.,  5., 10.,  1.33333333])
     2477        exact_edge_values = num.array([[  2.5,   2.,    1.5],
     2478                                       [  5.,   5.,   5., ],
     2479                                       [ 10.,   10.,  10. ],
     2480                                       [  3.,    2.,    0. ]])
     2481
     2482
     2483       
     2484        assert num.allclose(quantity.vertex_values,
     2485                            exact_vertex_values)
     2486        assert num.allclose(quantity.centroid_values, exact_centroid_values) #Centroid
     2487        assert num.allclose(quantity.edge_values, exact_edge_values)
     2488
     2489    def test_minimum(self):
     2490        quantity = Quantity(self.mesh4)
     2491
     2492        # get referece to data arrays
     2493        centroid_values = quantity.centroid_values
     2494        vertex_values = quantity.vertex_values
     2495        edge_values = quantity.edge_values
     2496
     2497        quantity.set_values([[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]],
     2498                            location = 'vertices')
     2499        assert num.allclose(quantity.vertex_values,
     2500                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
     2501
     2502        assert id(vertex_values) == id(quantity.vertex_values)
     2503        assert num.allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
     2504        assert num.allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
     2505                                                   [5., 5., 5.],
     2506                                                   [4.5, 4.5, 0.],
     2507                                                   [3.0, -1.5, -1.5]])
     2508
     2509
     2510
     2511        other_quantity = Quantity(self.mesh4)
     2512
     2513        other_quantity.set_values([[0,0,0], [1,1,6], [10,10,10], [0,0,4]],
     2514                            location = 'vertices')
     2515
     2516
     2517        #===============================
     2518        quantity.minimum(other_quantity)
     2519        #===============================
     2520
     2521        exact_vertex_values = num.array([[  0.,   0.,   0.],
     2522                                         [  1.,   1.,  5.],
     2523                                         [  0.,   0.,  9.],
     2524                                         [ -6.,   0.,  3.]])
     2525        exact_centroid_values = num.array([  0.,  2.66666667, 3.,  0.])
     2526        exact_edge_values = num.array([[  0.,   0.,   0.],
     2527                                       [  3.5,   3.5,   1., ],
     2528                                       [  4.5,   4.5,   0. ],
     2529                                       [  2.,   -1.5,  -1.5 ]])
     2530
     2531
     2532        assert num.allclose(quantity.vertex_values,
     2533                            exact_vertex_values)
     2534        assert num.allclose(quantity.centroid_values, exact_centroid_values) #Centroid
     2535        assert num.allclose(quantity.edge_values, exact_edge_values)
    24402536
    24412537
  • trunk/anuga_core/source/anuga/operators/set_stage_operator.py

    r8928 r8930  
    11"""
    2 Set value operators
     2Set stage operator
    33
    44Constraints: See GPL license in the user guide
     
    1717from anuga.geometry.polygon import inside_polygon
    1818
    19 from anuga.operators.base_operator import Operator
    20 from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    21                                               Modeltime_too_late
     19from anuga.operators.set_quantity_operator import Set_quantity_operator
     20
     21
    2222from anuga import indent
    2323
    2424
    2525
    26 class Set_stage_operator(Operator):
     26class Set_stage_operator(Set_quantity_operator):
    2727    """
    2828    Set the stage in a region
     
    3434    """
    3535
     36    get_stage = Set_quantity_operator.get_value
     37    set_stage = Set_quantity_operator.set_value
     38
    3639    def __init__(self,
    3740                 domain,
    3841                 stage=None,
    3942                 indices=None,
     43                 polygon=None,
     44                 center=None,
     45                 radius=None,
    4046                 description = None,
    4147                 label = None,
     
    4450
    4551
    46         Operator.__init__(self, domain, description, label, logging, verbose)
    47 
    48         #------------------------------------------
    49         # Local variables
    50         #------------------------------------------
    51         self.stage = stage
    52         self.indices = indices
    53 
    54 
    55     def __call__(self):
    56         """
    57         Apply rate to those triangles defined in indices
    58 
    59         indices == [], then don't apply anywhere
    60         indices == None, then apply everywhere
    61         otherwise apply for the specific indices
    62         """
    63 
    64 
    65         if self.indices is []:
    66             return
    67 
    68         stage = self.get_stage()
    69 
    70         if stage is None:
    71             return
    72 
    73         if self.verbose is True:
    74             log.critical('Stage of %s at time = %.2f = %f'
    75                          % (self.quantity_name, domain.get_time(), stage))
    76 
    77         if self.indices is None:
    78             self.stage_c[:] = stage
    79         else:
    80             self.stage_c[self.indices] = stage
    81 
    82 
    83     def get_stage(self, t=None):
    84         """Get value of stage at time t.
    85         If t not specified, return stage at current domain time
    86         """
    87 
    88         if t is None:
    89             t = self.domain.get_time()
    90 
    91         if callable(self.stage):
    92             try:
    93                 stage = self.stage(t)
    94             except Modeltime_too_early, e:
    95                 raise Modeltime_too_early(e)
    96             except Modeltime_too_late, e:
    97                 msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
    98                 msg += 'You can specify keyword argument default_rate in the '
    99                 msg += 'stage function to tell it what to do in the absence of time data.'
    100                 raise Modeltime_too_late(msg)
    101         else:
    102             stage = self.stage
    103 
    104 
    105         if stage is None:
    106             msg = ('Attribute stage must be specified in '+self.__name__+
    107                    ' before attempting to call it')
    108             raise Exception(msg)
    109 
    110         return stage
     52        Set_quantity_operator.__init__(self, domain,
     53                                       quantity = 'stage',
     54                                       value = stage,
     55                                       indices = indices,
     56                                       polygon = polygon,
     57                                       center = center,
     58                                       radius = radius,
     59                                       description = description,
     60                                       label = label,
     61                                       logging = logging,
     62                                       verbose = verbose)
    11163
    11264
     
    152104        assert radius is not None
    153105
    154 
    155         # Determine indices in update region
    156         N = domain.get_number_of_triangles()
    157         points = domain.get_centroid_coordinates(absolute=True)
    158 
    159 
    160         indices = []
    161 
    162         c = center
    163         r = radius
    164 
    165         self.center = center
    166         self.radius = radius
    167 
    168         intersect = False
    169         for k in range(N):
    170             x, y = points[k,:]    # Centroid
    171 
    172             if ((x-c[0])**2+(y-c[1])**2) < r**2:
    173                 intersect = True
    174                 indices.append(k)
    175 
    176 
    177         msg = 'No centroids intersect circle center'+str(center)+' radius '+str(radius)
    178         assert intersect, msg
    179 
    180 
    181 
    182 
    183         # It is possible that circle doesn't intersect with mesh (as can happen
    184         # for parallel runs
    185 
    186 
    187106        Set_stage_operator.__init__(self,
    188107                                    domain,
    189108                                    stage=stage,
    190                                     indices=indices,
     109                                    center=center,
     110                                    radius=radius,
    191111                                    verbose=verbose)
    192 
    193 
    194 
    195 
    196112
    197113#===============================================================================
     
    213129
    214130
    215         # Determine indices in update region
    216         N = domain.get_number_of_triangles()
    217         points = domain.get_centroid_coordinates(absolute=True)
    218 
    219 
    220         indices = inside_polygon(points, polygon)
    221         self.polygon = polygon
    222 
    223         # It is possible that circle doesn't intersect with mesh (as can happen
    224         # for parallel runs
    225 
    226131
    227132        Set_stage_operator.__init__(self,
    228133                               domain,
    229134                               stage=stage,
    230                                indices=indices,
     135                               polygon=polygon,
    231136                               verbose=verbose)
    232137
  • trunk/anuga_core/source/anuga/operators/test_set_elevation_operator.py

    r8928 r8930  
    869869           
    870870if __name__ == "__main__":
    871     suite = unittest.makeSuite(Test_set_elevation_operators, 'test')
     871    suite = unittest.makeSuite(Test_set_elevation_operator, 'test')
    872872    runner = unittest.TextTestRunner(verbosity=1)
    873873    runner.run(suite)
  • trunk/anuga_core/source/anuga/operators/test_set_stage_operator.py

    r8928 r8930  
    1313from anuga.config import time_format
    1414
    15 from set_stage_operators import *
     15from set_stage_operator import *
    1616
    1717import numpy as num
     
    253253
    254254        polygon = [(0.5,0.5), (1.5,0.5), (1.5,1.5), (0.5,1.5)]
     255        #operator = Polygonal_set_stage_operator(domain, stage=stage, polygon=polygon)
    255256        operator = Polygonal_set_stage_operator(domain, stage=stage, polygon=polygon)
    256257
  • trunk/anuga_core/source/anuga/utilities/function_utils.py

    r8920 r8930  
    2020    # or csllable
    2121    #------------------------------------------
    22     msg = "Input argument must be a scalar, or a function"
     22    msg = "Input argument must be a scalar, or a function or None"
    2323
     24    if function is None:
     25        return None
     26   
    2427    assert (isinstance(function, (int, float)) or
    2528            callable(function)), msg
  • trunk/anuga_core/source/anuga_validation_tests/case_studies/patong/README.txt

    r8908 r8930  
    2828-- Debugged it, and got rid of 'setup_model.py' to simplify the code
    2929
    30 -- God it running in parallel
     30-- Got it running in parallel
Note: See TracChangeset for help on using the changeset viewer.