Changeset 8483


Ignore:
Timestamp:
Jul 30, 2012, 5:38:32 PM (13 years ago)
Author:
steve
Message:

Very basic erosion operator

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

Legend:

Unmodified
Added
Removed

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

    r8482 r8483  
    1717
    1818from anuga.operators.base_operator import Operator
    19 from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    20                                               Modeltime_too_late
     19
    2120from anuga import indent
    2221
    23 def lineno():
    24     """Returns the current line number in our program."""
    25     import inspect
    26     return inspect.currentframe().f_back.f_back.f_lineno
    27 
    28 
    29 
    30 def stage_elev_info(self):
    31     print 80*"="
    32 
    33     print 'In Evolve: line number ', lineno()
    34     import inspect
    35     print inspect.getfile(lineno)
    36 
    37     print 80*"="
    38     ind = num.array([ 976,  977,  978,  979,  980,  981,  982,  983, 1016, 1017, 1018,
    39              1019, 1020, 1021, 1022, 1023])
    40     elev_v = self.get_quantity('elevation').vertex_values
    41     stage_v = self.get_quantity('stage').vertex_values
    42     elev_c = self.get_quantity('elevation').centroid_values
    43     stage_c = self.get_quantity('stage').centroid_values
    44 
    45     from pprint import pprint
    46     print 'elev_v, elev_c, elev_avg \n'
    47     pprint( num.concatenate( (elev_v[ind], (elev_c[ind]).reshape(16,1),
    48                                num.mean(elev_v[ind],axis=1).reshape(16,1)), axis = 1))
    49     print 'stage_v, stage_c, stage_avg \n'
    50     pprint( num.concatenate( (stage_v[ind], (stage_c[ind]).reshape(16,1),
    51                                num.mean(stage_v[ind],axis=1).reshape(16,1)), axis = 1))
    52 
    53    
    54 
    55     print 80*"="
    5622
    5723
    … …  
    8450        self.indices = indices
    8551
     52
     53
     54        self.elev_v  = self.domain.quantities['elevation'].vertex_values
     55        self.elev_e = self.domain.quantities['elevation'].edge_values
     56        self.stage_v = self.domain.quantities['stage'].vertex_values
     57
    8658        #------------------------------------------
    8759        # Need to turn off this optimization as it
    … …  
    9163        self.domain.optimise_dry_cells = 0
    9264
    93         #------------------------------------------
    94         # Find vertex nodes
    95         #------------------------------------------
    96         node_ids = set()
    97 
    98         for ind in self.indices:
    99             for k in [0,1,2]:
    100                 node_ids.add(self.domain.triangles[ind,k])
    101 
    102         self.node_ids = [ id for id in node_ids ]
    103 
    104 
    105         node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
    106 
    107         k = 0
    108         node_index[0] = 0
    109         for i in range(self.domain.number_of_nodes):
    110             node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
    111 
    112         self.node_index = node_index
    113 
    114         vertex_ids =[]
    115         for nid in self.node_ids:
    116             print nid,self.domain.number_of_triangles_per_node[nid]
    117             for vid in range(node_index[nid], node_index[nid+1]):
    118                 vidd = self.domain.vertex_value_indices[vid]
    119                 vertex_ids.append(vidd)
    120                 print '   ',nid, vid, vidd, vidd/3, vidd%3
    121 
    122         self.vol_ids  = num.array(vertex_ids,dtype=num.int)/3
    123         self.vols = num.array(list(set(self.vol_ids)), dtype=num.int)
    124         self.vert_ids = num.array(vertex_ids,dtype=num.int)%3
    125 
    126         print 'noe', self.domain.number_of_elements
    127         print 'non', self.domain.number_of_nodes
    128         #print self.vols
    129         #print self.domain.vertex_value_indic
    130         #print self.domain.number_of_triangles_per_node
    131         #print self.node_index
    132 
    133         print 'self.node_ids'
    134         print self.node_ids
    135 
    136         print 'self.indices'
    137         print self.indices
    138         #print self.domain.triangles[self.indices]
    139         #print self.vertex_ids
    140 
    141         self.dump_triangulation()
     65        #-----------------------------------------
     66        # Extra structures to support maintaining
     67        # continuity of elevation
     68        #-----------------------------------------
     69        self.setup_node_structures()
     70       
     71
     72    def update_quantities(self):
     73        """Update the vertex values of the quantities to model erosion
     74        """
     75
     76        t = self.get_time()
     77        dt = self.get_timestep()
     78
     79
     80        updated = True
     81
     82        if self.indices is None:
     83
     84            #--------------------------------------
     85            # Update all three vertices for each cell
     86            #--------------------------------------
     87            self.elev_v[:] = self.elev_v + 0.0
     88
     89        else:
     90
     91            #--------------------------------------
     92            # Update all three vertices for each cell
     93            #--------------------------------------
     94            ind = self.indices
     95            m = num.sqrt(self.xmom_c[ind]**2 + self.ymom_c[ind]**2)
     96            m = num.vstack((m,m,m)).T
     97            m = num.where(m>self.threshold, m, 0.0)
     98            self.elev_v[ind] = num.maximum(self.elev_v[ind] - m*dt, 0.0)
     99             #num.maximum(self.elev_v[ind] - momentum*dt, Z)
     100
     101
     102        return updated
     103
     104
     105
    142106
    143107    def __call__(self):
    … …  
    151115
    152116
    153 
    154117        if self.indices is []:
    155118            return
    156119
    157 
    158 
    159         #elevation = self.get_elevation()
    160 
    161 #        if self.verbose is True:
    162 #            log.critical('Bed of %s at time = %.2f = %f'
    163 #                         % (self.quantity_name, domain.get_time(), elevation))
    164 
    165         #if self.indices is None:
    166         #    self.elev_c[:] = elevation
    167         #else:
    168         #    self.elev_c[self.indices] = elevation
    169 
    170         t = self.get_time()
    171         dt = self.get_timestep()
    172 
    173 
    174         #print self.indices
    175        
    176 
    177         self.elev_v  = self.domain.quantities['elevation'].vertex_values
    178         self.stage_v = self.domain.quantities['stage'].vertex_values
    179 
    180 
    181         # Need to store water heights before change to ensure
    182         # no water lost or produced
    183         height_c = self.stage_c - self.elev_c
    184 
    185 
    186         #stage_elev_info(self.domain)
    187         #--------------------------------------------
    188         # Here we do the actual erosion
    189         #--------------------------------------------
    190 
    191         if t < 10.0 and t > 7.0:
    192             if self.indices is None:
    193                 self.elev_v[:] = self.elev_v + 0.0
    194             else:
    195                 self.elev_v[self.indices] -= 0.1*dt
    196 
    197 
    198                 #self.elev_v[self.vol_ids, self.vert_ids] = \
    199                 #   num.maximum(self.elev_v[self.vol_ids, self.vert_ids] - 0.1*dt, 0.0)
    200 
    201                 #self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1)
    202 
    203         #for nid in self.node_id:
    204         #    print 'nid ', nid
    205         #    print 'vvi ', self.domain.vertex_value_indices[nid]
    206 
    207 
    208         #stage_elev_info(self.domain)
    209        
    210 
    211         # FIXME SR: At present need to ensure the elevation is continuous
    212         # In future with discontinuous bed we will not need to do this.
    213         self.domain.quantities['elevation'].smooth_vertex_values()
    214 
    215         # Need to do this faster.
    216 
    217         #stage_elev_info(self.domain)
    218 
    219         self.domain.quantities['elevation'].interpolate()
    220 
    221         #stage_elev_info(self.domain)
    222 
    223         #self.elev_c = self.domain.quantities['elevation'].centroid_values
    224 
    225 #        # Fix up water conservation
    226         self.stage_c[:] = self.elev_c +  height_c
    227 
    228         #stage_elev_info(self.domain)
    229 
    230 
    231         #old_flag = self.domain.optimise_dry_cells
    232         #self.domain.optimise_dry_cells = 0
    233         #self.domain.distribute_to_vertices_and_edges()
    234         #self.domain.optimise_dry_cells = old_flag
    235 
    236         #stage_elev_info(self.domain)
    237         #print 'time in erosion ',self.get_time(), dt
     120        #------------------------------------------
     121        # Apply changes to elevation vertex values
     122        # via the update_quantites routine
     123        #------------------------------------------
     124        if not self.update_quantities():
     125            return
     126
     127
     128        #------------------------------------------
     129        # Cleanup elevation and stage quantity values
     130        #-----------------------------------------
     131        if self.indices is None:
     132
     133            #--------------------------------------
     134            # Update all three vertices for each cell
     135            #--------------------------------------
     136            self.elev_v[:] = self.elev_v + 0.0
     137
     138            #--------------------------------------
     139            # Make elevation continuous and clean up
     140            # stage values to ensure conservation
     141            #--------------------------------------
     142            height_c = self.stage_c - self.elev_c
     143            self.domain.quantities['elevation'].smooth_vertex_values()
     144            self.domain.quantities['elevation'].interpolate()
     145            self.stage_c[:] = self.elev_c +  height_c
     146
     147
     148        else:
     149
     150            #--------------------------------------
     151            # Make elevation continuous and clean up
     152            # stage values to ensure conservation
     153            #--------------------------------------
     154            height_c = self.stage_c[self.vols] - self.elev_c[self.vols]
     155            for nid in self.node_ids:
     156                non = self.domain.number_of_triangles_per_node[nid]
     157
     158                vid = num.arange(self.node_index[nid], self.node_index[nid+1],dtype=num.int)
     159                vidd = self.domain.vertex_value_indices[vid]
     160
     161                self.elev_v[vidd/3,vidd%3] = num.sum(self.elev_v[vidd/3,vidd%3])/non
     162
     163
     164            #--------------------------------------
     165            # clean up the centroid values and edge values
     166            #--------------------------------------
     167            self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1)
     168
     169            self.elev_e[self.vols,0] = 0.5*(self.elev_v[self.vols,1]+ self.elev_v[self.vols,2])
     170            self.elev_e[self.vols,1] = 0.5*(self.elev_v[self.vols,2]+ self.elev_v[self.vols,0])
     171            self.elev_e[self.vols,2] = 0.5*(self.elev_v[self.vols,0]+ self.elev_v[self.vols,1])
     172
     173            self.stage_c[self.vols] = self.elev_c[self.vols] +  height_c
     174
    238175
    239176
    … …  
    263200        try:
    264201            import matplotlib
    265             matplotlib.use('Agg')
     202            #matplotlib.use('Agg')
    266203            import matplotlib.pyplot as plt
    267204            import matplotlib.tri as tri
    … …  
    312249        noe = self.domain.number_of_elements
    313250
     251
    314252        fx = vertices[:,0].reshape(noe,3)
    315253        fy = vertices[:,1].reshape(noe,3)
    316254
    317255
    318         fx = fx[self.indices].flatten()
    319         fy = fy[self.indices].flatten()
    320 
    321         print 'fx', fx.shape
     256
     257
     258        #-------------------------------------------
     259        # Neighbourhood Area
     260        #-------------------------------------------
     261        fx1 = fx[self.vols].flatten()
     262        fy1 = fy[self.vols].flatten()
     263
     264        print 'fx1', fx1.shape
     265
     266        print self.vols
     267        #gx = vertices[ghost_mask,0]
     268        #gy = vertices[ghost_mask,1]
     269
     270
     271        ## Plot full triangles
     272        n = int(len(fx1)/3)
     273        triang = num.array(range(0,3*n))
     274        triang.shape = (n, 3)
     275        print triang
     276        plt.triplot(fx1, fy1, triang, 'go-')
     277
     278
     279
     280        self.vols
     281        #plt.plot()
     282
     283        #-------------------------------------------
     284        # Update Area
     285        #-------------------------------------------
     286        fx0 = fx[self.indices].flatten()
     287        fy0 = fy[self.indices].flatten()
     288
     289        print 'fx0', fx0.shape
    322290
    323291        print self.indices
    … …  
    327295
    328296        ## Plot full triangles
    329         n = int(len(fx)/3)
    330 
    331         Z = num.ones((3*n,),dtype=num.float)
    332         print Z.shape
    333 
     297        n = int(len(fx0)/3)
    334298        triang = num.array(range(0,3*n))
    335299        triang.shape = (n, 3)
    336300        print triang
    337         plt.triplot(fx, fy, triang, 'o-')
    338         plt.tripcolor(fx,fy, triang, Z)
     301        plt.triplot(fx0, fy0, triang, 'bo-')
     302
     303
     304        #-------------------------------------------
     305        # Update Nodes
     306        #-------------------------------------------
     307        fx0 = fx[self.indices].flatten()
     308        fy0 = fy[self.indices].flatten()
     309
     310        print 'fx0', fx0.shape
     311
     312        print self.indices
     313        #gx = vertices[ghost_mask,0]
     314        #gy = vertices[ghost_mask,1]
     315
     316
     317        ## Plot full triangles
     318        n = int(len(fx0)/3)
     319        triang = num.array(range(0,3*n))
     320        triang.shape = (n, 3)
     321        print triang
     322        plt.triplot(fx0, fy0, triang, 'bo-')
     323
     324
     325
     326
     327        fx2 = fx[self.vol_ids,self.vert_ids]
     328        fy2 = fy[self.vol_ids,self.vert_ids]
     329
     330        print 'fx2', fx2.shape
     331
     332        plt.plot(fx2,fy2,'yo')
     333
     334        #plt.tripcolor(fx,fy, triang, Z)
    339335       
    340336        ## Plot ghost triangles
    … …  
    352348
    353349
     350    def setup_node_structures(self):
     351        """ For setting elevation we need to
     352        ensure that the elevation quantity remains
     353        continuous (at least for version 1.3 of anuga)
     354
     355        So we need to find all the vertices that need to
     356        update within each timestep.
     357        """
     358
     359        node_ids = set()
     360
     361        for ind in self.indices:
     362            for k in [0,1,2]:
     363                node_ids.add(self.domain.triangles[ind,k])
     364
     365        self.node_ids = [ id for id in node_ids ]
     366
     367
     368        node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
     369
     370        k = 0
     371        node_index[0] = 0
     372        for i in range(self.domain.number_of_nodes):
     373            node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
     374
     375        self.node_index = node_index
     376
     377        vertex_ids =[]
     378        for nid in self.node_ids:
     379            #print nid,self.domain.number_of_triangles_per_node[nid]
     380            for vid in range(node_index[nid], node_index[nid+1]):
     381                vidd = self.domain.vertex_value_indices[vid]
     382                vertex_ids.append(vidd)
     383                #print '   ',nid, vid, vidd, vidd/3, vidd%3
     384
     385        self.vol_ids  = num.array(vertex_ids,dtype=num.int)/3
     386        self.vols = num.array(list(set(self.vol_ids)), dtype=num.int)
     387        self.vert_ids = num.array(vertex_ids,dtype=num.int)%3
     388
    354389
    355390
    … …  
    453488
    454489
     490def lineno():
     491    """Returns the current line number in our program."""
     492    import inspect
     493    return inspect.currentframe().f_back.f_back.f_lineno
     494
     495
     496
     497def stage_elev_info(self):
     498    print 80*"="
     499
     500    print 'In Evolve: line number ', lineno()
     501    import inspect
     502    print inspect.getfile(lineno)
     503
     504    print 80*"="
     505    ind = num.array([ 976,  977,  978,  979,  980,  981,  982,  983, 1016, 1017, 1018,
     506             1019, 1020, 1021, 1022, 1023])
     507    elev_v = self.get_quantity('elevation').vertex_values
     508    stage_v = self.get_quantity('stage').vertex_values
     509    elev_c = self.get_quantity('elevation').centroid_values
     510    stage_c = self.get_quantity('stage').centroid_values
     511
     512    from pprint import pprint
     513    print 'elev_v, elev_c, elev_avg \n'
     514    pprint( num.concatenate( (elev_v[ind], (elev_c[ind]).reshape(16,1),
     515                               num.mean(elev_v[ind],axis=1).reshape(16,1)), axis = 1))
     516    print 'stage_v, stage_c, stage_avg \n'
     517    pprint( num.concatenate( (stage_v[ind], (stage_c[ind]).reshape(16,1),
     518                               num.mean(stage_v[ind],axis=1).reshape(16,1)), axis = 1))
     519
     520    print 80*"="

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

    r8482 r8483  
    2020length = 24.
    2121width = 5.
    22 dx = dy = 1.0 #.1           # Resolution: Length of subdivisions on both axes
     22dx = dy = 0.1 #.1           # Resolution: Length of subdivisions on both axes
    2323
    2424points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
    2525                                               len1=length, len2=width)
    2626domain = Domain(points, vertices, boundary)
    27 domain.set_name('erosion_polygon') # Output name
     27domain.set_name('polygon_erosion') # Output name
    2828print domain.statistics()
    2929domain.set_quantities_to_be_stored({'elevation': 2,
    … …  
    8484#------------------------------------------------------------------------------
    8585
    86 for t in domain.evolve(yieldstep=0.2, finaltime=2.0):
     86for t in domain.evolve(yieldstep=0.2, finaltime=30.0):
    8787    domain.print_timestepping_statistics()
    8888    domain.print_operator_timestepping_statistics()

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

    r8477 r8483  
    8282        return 0.0
    8383    else:
    84         return 0.05
     84        return 1.0
    8585
    8686
    … …  
    114114    domain.print_operator_timestepping_statistics()
    115115
    116     #w = domain.get_quantity('stage').\
    117     #    get_values(interpolation_points=[[18, 2.5]])
    118     #print 'Level at gauge point = %.2fm' % w
    119 
    120     #z = domain.get_quantity('elevation').\
    121     #    get_values(interpolation_points=[[12, 3]])
    122     #print 'Elevation at pole location = %.2fm' % z
    123 
    124 #    # Start variable elevation after 10 seconds
    125 #    if t > 10 and not (shrinking or growing or done):
    126 #        growing = True
    127 #
    128 #    # Stop growing when pole has reached a certain height
    129 #    if t > 16 and growing:
    130 #        growing = False
    131 #        shrinking = False
    132 #
    133 #    # Start shrinking
    134 #    if t > 20:
    135 #        shrinking = True
    136 #        growing = False
    137 #
    138 #    # Stop changing when pole has shrunk to original level
    139 #    if t > 25 and shrinking:
    140 #        done = True
    141 #        shrinking = growing = False
    142 #        domain.set_quantity('elevation', topography)
    143 #
    144 #    # Grow or shrink
    145 #    if growing:
    146 #        domain.add_quantity('elevation', pole_increment)
    147 #
    148 #    if shrinking:
    149 #        domain.add_quantity('elevation', lambda x,y: -2*pole_increment(x,y))
    150116
    151117
    152118
    153 
    154 

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

    r8477 r8483  
    1717
    1818from anuga.operators.base_operator import Operator
    19 from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
    20                                               Modeltime_too_late
     19
    2120from anuga import indent
    2221
    … …  
    5049        self.elevation = elevation
    5150        self.indices = indices
     51       
     52        #------------------------------------------
     53        # Extra aliases for changing elevation at
     54        # vertices and edges
     55        #------------------------------------------
     56        self.elev_v  = self.domain.quantities['elevation'].vertex_values
     57        self.elev_e = self.domain.quantities['elevation'].edge_values
     58
     59        #------------------------------------------
     60        # Need to turn off this optimization as it
     61        # doesn't fixup the relationship between
     62        # bed and stage vertex values in dry region
     63        #------------------------------------------
     64        self.domain.optimise_dry_cells = 0
     65
     66        #-----------------------------------------
     67        # Extra structures to support maintaining
     68        # continuity of elevation
     69        #-----------------------------------------
     70        self.setup_node_structures()
    5271
    5372
    … …  
    6180        """
    6281
    63         #if self.indices is []:
    64         #    return
    65 
    66         #elevation = self.get_elevation()
    67 
    68 #        if self.verbose is True:
    69 #            log.critical('Bed of %s at time = %.2f = %f'
    70 #                         % (self.quantity_name, domain.get_time(), elevation))
    71 
    72         #if self.indices is None:
    73         #    self.elev_c[:] = elevation
    74         #else:
    75         #    self.elev_c[self.indices] = elevation
    76 
    77         t = self.get_time()
    78         dt = self.get_timestep()
    79 
    80         v_coors = self.domain.vertex_coordinates
    81         self.elev_v = self.domain.quantities['elevation'].vertex_values
    82 
    83 
     82        if self.indices is []:
     83            return
     84
     85        #------------------------------------------
     86        # Apply changes to elevation vertex values
     87        # via the update_quantites routine
     88        #------------------------------------------
     89        if not self.update_quantities():
     90            return
     91
     92
     93        #------------------------------------------
     94        # Cleanup elevation and stage quantity values
     95        #-----------------------------------------
    8496        if self.indices is None:
    85             self.elev_v[:] = self.elev_v + 0.0
     97
     98            #--------------------------------------
     99            # Make elevation continuous and clean up
     100            # stage values to ensure conservation
     101            #--------------------------------------
     102            height_c = self.stage_c - self.elev_c
     103            self.domain.quantities['elevation'].smooth_vertex_values()
     104            self.domain.quantities['elevation'].interpolate()
     105            self.stage_c[:] = self.elev_c +  height_c
     106
     107
    86108        else:
    87             self.elev_v[self.indices] += self.elevation(t)*dt
    88 
    89         ### make sure centroid is correct as well
    90        
    91         #self.domain.add_quantity('elevation', lambda x,y: dt*self.elevation(x,y,t))
    92 
    93 
    94 
    95         # clean up discontinuities for now
    96         #self.domain.quantities['elevation'].smooth_vertex_values()
    97 
    98 
    99 
     109
     110            #--------------------------------------
     111            # Make elevation continuous and clean up
     112            # stage values to ensure conservation
     113            #--------------------------------------
     114            height_c = self.stage_c[self.vols] - self.elev_c[self.vols]
     115            for nid in self.node_ids:
     116                non = self.domain.number_of_triangles_per_node[nid]
     117
     118                vid = num.arange(self.node_index[nid], self.node_index[nid+1],dtype=num.int)
     119                vidd = self.domain.vertex_value_indices[vid]
     120
     121                self.elev_v[vidd/3,vidd%3] = num.sum(self.elev_v[vidd/3,vidd%3])/non
     122
     123
     124            #--------------------------------------
     125            # clean up the centroid values and edge values
     126            #--------------------------------------
     127            self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1)
     128
     129            self.elev_e[self.vols,0] = 0.5*(self.elev_v[self.vols,1]+ self.elev_v[self.vols,2])
     130            self.elev_e[self.vols,1] = 0.5*(self.elev_v[self.vols,2]+ self.elev_v[self.vols,0])
     131            self.elev_e[self.vols,2] = 0.5*(self.elev_v[self.vols,0]+ self.elev_v[self.vols,1])
     132
     133            self.stage_c[self.vols] = self.elev_c[self.vols] +  height_c
     134
     135
     136
     137    def update_quantities(self):
     138        """Update the vertex values of the quantities to model erosion
     139        """
     140
     141
     142        elevation = self.get_elevation()
     143
     144        updated = True
     145
     146        if self.indices is None:
     147
     148            #--------------------------------------
     149            # Update all three vertices for each cell
     150            #--------------------------------------
     151            self.elev_v[:] = elevation
     152
     153        else:
     154
     155            #--------------------------------------
     156            # Update all three vertices for each cell
     157            #--------------------------------------
     158            self.elev_v[self.indices] = elevation
     159
     160
     161        return updated
    100162
    101163    def get_elevation(self, t=None):
    … …  
    104166        """
    105167
     168        from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
     169                                                      Modeltime_too_late
     170                                                     
    106171        if t is None:
    107172            t = self.domain.get_time()
    … …  
    128193        return elevation
    129194
     195    def setup_node_structures(self):
     196        """ For setting elevation we need to
     197        ensure that the elevation quantity remains
     198        continuous (at least for version 1.3 of anuga)
     199
     200        So we need to find all the vertices that need to
     201        update within each timestep.
     202        """
     203
     204        node_ids = set()
     205
     206        for ind in self.indices:
     207            for k in [0,1,2]:
     208                node_ids.add(self.domain.triangles[ind,k])
     209
     210        self.node_ids = [ id for id in node_ids ]
     211
     212
     213        node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
     214
     215        k = 0
     216        node_index[0] = 0
     217        for i in range(self.domain.number_of_nodes):
     218            node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
     219
     220        self.node_index = node_index
     221
     222        vertex_ids =[]
     223        for nid in self.node_ids:
     224            #print nid,self.domain.number_of_triangles_per_node[nid]
     225            for vid in range(node_index[nid], node_index[nid+1]):
     226                vidd = self.domain.vertex_value_indices[vid]
     227                vertex_ids.append(vidd)
     228                #print '   ',nid, vid, vidd, vidd/3, vidd%3
     229
     230        self.vol_ids  = num.array(vertex_ids,dtype=num.int)/3
     231        self.vols = num.array(list(set(self.vol_ids)), dtype=num.int)
     232        self.vert_ids = num.array(vertex_ids,dtype=num.int)%3
     233
    130234
    131235
Note: See TracChangeset for help on using the changeset viewer.