Changeset 8928


Ignore:
Timestamp:
Jun 25, 2013, 7:34:21 PM (12 years ago)
Author:
steve
Message:

Small change to set operators

Location:
trunk/anuga_core/source/anuga/operators
Files:
1 edited
4 moved

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/operators/run_set_elevation.py

    r8873 r8928  
    3636    """Complex topography defined by a function of vectors x and y."""
    3737
     38
    3839    z = -x/100
    3940
    40     N = len(x)
    41     for i in range(N):
    42         # Step
    43         if 2 < x[i] < 4:
    44             z[i] += 0.4 - 0.05*y[i]
    4541
    46         # Permanent pole
    47         if (x[i] - 8)**2 + (y[i] - 2)**2 < 0.4**2:
    48             z[i] += 1
     42    # Step
     43    id_s = (2 < x) & (x < 4)
     44
     45    z[id_s] += 0.4 - 0.05*y[id_s]
     46
     47    # Permanent pole
     48    id_p = (x - 8)**2 + (y - 2)**2 < 0.4**2
     49
     50    z[id_p] += 1
     51
     52           
    4953    return z
    5054
    … …  
    8185    if t<10:
    8286        return 0.0
    83     elif t>15:
     87    elif t>12:
    8488        return 0.0
    8589    else:
    … …  
    105109
    106110from anuga.operators.set_elevation_operators import Set_elevation_operator
    107 
    108111op1 = Set_elevation_operator(domain, elevation=pole, radius=0.5, center = (12.0,3.0))
    109112
    110 from anuga.operators.set_elevation_operators import Set_elevation
    111 op2 = Set_elevation(domain, elevation = topography, radius = 1.0, center = (12.0,3.0))
     113from anuga.operators.set_elevation import Set_elevation
     114op2 = Set_elevation(domain, elevation = topography)
    112115
    113116growing = False
    … …  
    119122    domain.print_operator_timestepping_statistics()
    120123
    121     if num.allclose(17.0, t):
     124    if num.allclose(15.0, t):
     125        op1.set_elevation(None)
    122126        op2()
    123127

  • trunk/anuga_core/source/anuga/operators/set_elevation_oprator.py

    r8926 r8928  
    1717
    1818from anuga.operators.base_operator import Operator
     19from anuga.operators.set_elevation import Set_elevation
    1920
    2021from anuga import indent
    21 
    22 class Set_elevation(object):
    23     """
    24     Helper class to setup calculation of elevation
    25     associated with a region (defined by indices, polygon or center/radius
    26     """
    27 
    28     def __init__(self,
    29                  domain,
    30                  elevation=None,
    31                  indices=None,
    32                  polygon=None,
    33                  center=None,
    34                  radius=None):
    35 
    36         #------------------------------------------
    37         # Local variables
    38         #------------------------------------------
    39         self.elevation = elevation
    40         self.indices = indices
    41         self.center = center
    42         self.radius = radius
    43         self.polygon = polygon
    44 
    45    
    46         #------------------------------------------
    47         #  Useful aliases
    48         #------------------------------------------
    49         self.domain = domain
    50         self.stage_c = self.domain.quantities['stage'].centroid_values
    51         self.xmom_c  = self.domain.quantities['xmomentum'].centroid_values
    52         self.ymom_c  = self.domain.quantities['ymomentum'].centroid_values
    53         self.elev_c  = self.domain.quantities['elevation'].centroid_values
    54         self.coord_c = self.domain.centroid_coordinates
    55         self.areas = self.domain.areas
    56 
    57 
    58         #------------------------------------------
    59         # Extra aliases for changing elevation at
    60         # vertices and edges
    61         #------------------------------------------
    62         self.elev_v  = self.domain.quantities['elevation'].vertex_values
    63         self.elev_e = self.domain.quantities['elevation'].edge_values
    64 
    65         #------------------------------------------
    66         # Need to turn off this optimization as it
    67         # doesn't fixup the relationship between
    68         # bed and stage vertex values in dry region
    69         #------------------------------------------
    70         self.domain.optimise_dry_cells = 0
    71 
    72 
    73         #-------------------------------------------
    74         # Work out region:
    75         #-------------------------------------------
    76         if self.indices is not None:
    77             # This overrides polygon, center and radius
    78             pass
    79         elif (self.center is not None) and (self.radius is not None):
    80 
    81             assert self.indices is None
    82             assert self.polygon is None
    83 
    84             self.setup_indices_circle()
    85 
    86         elif (self.polygon is not None):
    87 
    88             assert self.indices is None
    89 
    90             self.setup_indices_polygon()
    91         else:
    92             assert self.indices is None or self.indices is []
    93 
    94 
    95         #-----------------------------------------
    96         # Extra structures to support maintaining
    97         # continuity of elevation
    98         #-----------------------------------------
    99         self.setup_node_structures()
    100 
    101 
    102     def __call__(self):
    103         """
    104         Apply rate to those triangles defined in indices
    105 
    106         indices == [], then don't apply anywhere
    107         indices == None, then apply everywhere
    108         otherwise apply for the specific indices
    109         """
    110 
    111         if self.indices is []:
    112             return
    113 
    114         #------------------------------------------
    115         # Apply changes to elevation vertex values
    116         # via the update_quantites routine
    117         #------------------------------------------
    118         if not self.update_quantities():
    119             return
    120 
    121 
    122         #------------------------------------------
    123         # Cleanup elevation and stage quantity values
    124         #-----------------------------------------
    125         if self.indices is None:
    126 
    127             #--------------------------------------
    128             # Make elevation continuous and clean up
    129             # stage values to ensure conservation
    130             #--------------------------------------
    131             height_c = self.stage_c - self.elev_c
    132             self.domain.quantities['elevation'].smooth_vertex_values()
    133             self.domain.quantities['elevation'].interpolate()
    134             self.stage_c[:] = self.elev_c +  height_c
    135 
    136 
    137         else:
    138 
    139             #--------------------------------------
    140             # Make elevation continuous and clean up
    141             # stage values to ensure conservation
    142             #--------------------------------------
    143             height_c = self.stage_c[self.vols] - self.elev_c[self.vols]
    144             for nid in self.node_ids:
    145                 non = self.domain.number_of_triangles_per_node[nid]
    146 
    147                 vid = num.arange(self.node_index[nid], self.node_index[nid+1],dtype=num.int)
    148                 vidd = self.domain.vertex_value_indices[vid]
    149 
    150                 self.elev_v[vidd/3,vidd%3] = num.sum(self.elev_v[vidd/3,vidd%3])/non
    151 
    152 
    153             #--------------------------------------
    154             # clean up the centroid values and edge values
    155             #--------------------------------------
    156             self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1)
    157 
    158             self.elev_e[self.vols,0] = 0.5*(self.elev_v[self.vols,1]+ self.elev_v[self.vols,2])
    159             self.elev_e[self.vols,1] = 0.5*(self.elev_v[self.vols,2]+ self.elev_v[self.vols,0])
    160             self.elev_e[self.vols,2] = 0.5*(self.elev_v[self.vols,0]+ self.elev_v[self.vols,1])
    161 
    162             self.stage_c[self.vols] = self.elev_c[self.vols] +  height_c
    163 
    164 
    165 
    166     def update_quantities(self):
    167         """Update the vertex values of the quantities to model erosion
    168         """
    169 
    170 
    171         elevation = self.get_elevation()
    172 
    173         if elevation is None:
    174             return False
    175 
    176         updated = True
    177 
    178         if self.indices is None:
    179 
    180             #--------------------------------------
    181             # Update all three vertices for each cell
    182             #--------------------------------------
    183             self.elev_v[:] = elevation
    184 
    185         else:
    186 
    187             #--------------------------------------
    188             # Update all three vertices for each cell
    189             #--------------------------------------
    190             self.elev_v[self.indices] = elevation
    191 
    192 
    193         return updated
    194 
    195     def get_elevation(self, t=None):
    196         """Get value of elevation at time t.
    197         If t not specified, return elevation at current domain time
    198         """
    199 
    200         from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    201                                                       Modeltime_too_late
    202                                                      
    203         if t is None:
    204             t = self.domain.get_time()
    205 
    206         if callable(self.elevation):
    207             try:
    208                 elevation = self.elevation(t)
    209             except Modeltime_too_early, e:
    210                 raise Modeltime_too_early(e)
    211             except Modeltime_too_late, e:
    212                 msg = '%s: ANUGA is trying to run longer than specified data.\n' %str(e)
    213                 msg += 'You can specify keyword argument default_rate in the '
    214                 msg += 'rate operator to tell it what to do in the absence of time data.'
    215                 raise Modeltime_too_late(msg)
    216         else:
    217             elevation = self.elevation
    218 
    219 
    220 #        if elevation is None:
    221 #            msg  = 'Attribute elevation must be specified'
    222 #            raise Exception(msg)
    223 
    224         return elevation
    225 
    226 
    227 
    228     def setup_indices_circle(self):
    229 
    230         # Determine indices in circular region
    231         N = self.domain.get_number_of_triangles()
    232         points = self.domain.get_centroid_coordinates(absolute=True)
    233 
    234         indices = []
    235 
    236         c = self.center
    237         r = self.radius
    238 
    239 
    240         intersect = False
    241         for k in range(N):
    242             x, y = points[k,:]    # Centroid
    243 
    244             if ((x-c[0])**2+(y-c[1])**2) < r**2:
    245                 intersect = True
    246                 indices.append(k)
    247 
    248         self.indices = indices
    249 
    250         msg = 'No centroids intersect circle center'+str(c)+' radius '+str(r)
    251         assert intersect, msg
    252 
    253 
    254     def setup_indices_polygon(self):
    255 
    256         # Determine indices for polygonal region
    257         points = self.domain.get_centroid_coordinates(absolute=True)
    258 
    259 
    260         self.indices = inside_polygon(points, self.polygon)
    261 
    262 
    263 
    264 
    265     def setup_node_structures(self):
    266         """ For setting elevation we need to
    267         ensure that the elevation quantity remains
    268         continuous (at least for version 1.3 of anuga)
    269 
    270         So we need to find all the vertices that need to
    271         update within each timestep.
    272         """
    273 
    274 
    275 
    276         if self.indices is None or self.indices is []:
    277             self.vol_ids  = None
    278             self.vols = None
    279             self.vert_ids = None
    280             return
    281 
    282 
    283        
    284         node_ids = set()
    285 
    286         for ind in self.indices:
    287             for k in [0,1,2]:
    288                 node_ids.add(self.domain.triangles[ind,k])
    289 
    290         self.node_ids = [ id for id in node_ids ]
    291 
    292 
    293 
    294         node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
    295 
    296         # FIXME: SR Don't we calculate this for the domain already!
    297         k = 0
    298         node_index[0] = 0
    299         for i in range(self.domain.number_of_nodes):
    300             node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
    301 
    302         self.node_index = node_index
    303 
    304         vertex_ids =[]
    305         for nid in self.node_ids:
    306             #print nid,self.domain.number_of_triangles_per_node[nid]
    307             for vid in range(node_index[nid], node_index[nid+1]):
    308                 vidd = self.domain.vertex_value_indices[vid]
    309                 vertex_ids.append(vidd)
    310                 #print '   ',nid, vid, vidd, vidd/3, vidd%3
    311 
    312         self.vol_ids  = num.array(vertex_ids,dtype=num.int)/3
    313         self.vols = num.array(list(set(self.vol_ids)), dtype=num.int)
    314         self.vert_ids = num.array(vertex_ids,dtype=num.int)%3
    315 
    316 
    317 
    318 
    31922
    32023
    … …  
    37679        #message  = indent + self.label + ': Set_elevation = ' + str('')
    37780        #message  += ' at center '+str(self.center)
    378         return 'test'
     81
     82        return 'Set_elevation_operator active at time '+str(self.domain.get_time())
    37983
    38084

  • trunk/anuga_core/source/anuga/operators/test_set_elevation_operator.py

    r8926 r8928  
    1313from anuga.config import time_format
    1414
    15 from set_elevation_operators import *
     15from set_elevation_operator import *
    1616
    1717import numpy as num
    … …  
    2121
    2222
    23 class Test_set_elevation_operators(unittest.TestCase):
     23class Test_set_elevation_operator(unittest.TestCase):
    2424    def setUp(self):
    2525        pass
Note: See TracChangeset for help on using the changeset viewer.