Changeset 7968


Ignore:
Timestamp:
Aug 23, 2010, 2:48:30 PM (9 years ago)
Author:
steve
Message:

Commiting changes of refactor of culvert operator

Location:
trunk/anuga_core/source/anuga
Files:
3 edited

Legend:

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

    r7711 r7968  
    884884                return i
    885885
    886         return
     886        msg = 'Point %s not found within a triangle' %str(point)
     887        raise Exception, msg
     888
    887889
    888890
  • trunk/anuga_core/source/anuga/structures/culvert_operator.py

    r7962 r7968  
    2020class Above_interval(Exception): pass
    2121
    22    
     22
     23
     24class Inlet:
     25    """Contains information associated with each inlet
     26    """
     27
     28    def __init__(self,domain,polygon,enquiry_point,inlet_vector):
     29
     30        self.domain = domain
     31        self.domain_bounding_polygon = self.domain.get_boundary_polygon()
     32        self.polygon = polygon
     33        self.enquiry_point = enquiry_point
     34        self.inlet_vector = inlet_vector
     35
     36        # FIXME (SR) Using get_triangle_containing_point which needs to be sped up
     37        self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_point)
     38
     39        self.compute_inlet_triangle_indices()
     40
     41
     42    def compute_inlet_averages(self):
     43
     44
     45
     46        self.cell_indices = self.exchange_triangle_indices[0]
     47        self.areas    = areas[cell_indices]
     48
     49
     50        # Inlet Averages
     51        self.heights            = stage[cell_indices]-elevation[cell_indices]
     52        self.total_water_volume = num.sum(self.heights*self.areas)
     53        self.average_height     = self.total_water_volume/self.total_inlet_area
     54
     55
     56
     57
     58    def compute_inlet_triangle_indices(self):
     59
     60        # Get boundary (in absolute coordinates)
     61        bounding_polygon = self.domain_bounding_polygon
     62        centroids = self.domain.get_centroid_coordinates(absolute=True)
     63
     64        inlet_polygon = self.polygon
     65
     66        # Check that polygon lies within the mesh.
     67        for point in inlet_polygon:
     68                msg = 'Point %s in polygon for forcing term' % str(point)
     69                msg += ' did not fall within the domain boundary.'
     70                assert is_inside_polygon(point, bounding_polygon), msg
     71
     72
     73        self.triangle_indices = inside_polygon(centroids, inlet_polygon)
     74
     75        if len(self.triangle_indices) == 0:
     76            region = 'Inlet polygon=%s' % (inlet_polygon)
     77            msg = 'No triangles have been identified in '
     78            msg += 'specified region: %s' % region
     79            raise Exception, msg
     80
     81        # Compute exchange area as the sum of areas of triangles identified
     82        # by polygon
     83        self.area = 0.0
     84        for j in self.triangle_indices:
     85            self.area += self.domain.areas[j]
     86
     87
     88        msg = 'Inlet exchange area has area = %f' % self.area
     89        assert self.area > 0.0
     90
     91
    2392
    2493class Generic_box_culvert:
     
    60129        self.filename = None
    61130       
    62         # Create the fundamental culvert polygons from polygon
     131        # Create the fundamental culvert polygons and create inlet objects
    63132        self.create_culvert_polygons()
    64         self.compute_enquiry_indices()
    65         self.check_culvert_inside_domain()
    66         self.compute_exchange_triangle_indices()
    67 
    68 
    69         #self.print_stats()
     133
     134        #FIXME (SR) Put this into a foe loop to deal with more inlets
     135        self.inlets = []
     136        polygon0 = self.exchange_polygons[0]
     137        enquiry_pt0 = self.enquiry_points[0]
     138        inlet0_vector = self.culvert_vector
     139
     140        self.inlets.append(Inlet(self.domain,polygon0,enquiry_pt0,inlet0_vector))
     141
     142        polygon1 = self.exchange_polygons[1]
     143        enquiry_pt1 = self.enquiry_points[1]
     144        inlet1_vector = - self.culvert_vector
     145
     146        self.inlets.append(Inlet(self.domain,polygon1,enquiry_pt1, inlet1_vector))
     147
     148
     149        # aliases to quantity centroid values and cell areas
     150        self.areas      = self.domain.areas
     151        self.stage      = self.domain.quantities['stage'].centroid_values
     152        self.elevation  = self.domain.quantities['elevation'].centroid_values
     153        self.xmom       = self.domain.quantities['xmomentum'].centroid_values
     154        self.ymom       = self.domain.quantities['ymomentum'].centroid_values
     155
     156 
     157        self.print_stats()
    70158
    71159
     
    79167
    80168
    81         inlet_indices  = self.exchange_triangle_indices[0]
    82         outlet_indices = self.exchange_triangle_indices[1]
    83 
    84         areas      = self.domain.areas
    85         stage     = self.domain.quantities['stage'].centroid_values
    86         elevation = self.domain.quantities['elevation'].centroid_values
    87        
    88         xmom = self.domain.quantities['xmomentum'].centroid_values
    89         ymom = self.domain.quantities['ymomentum'].centroid_values
    90 
    91         # Inlet averages
    92         inlet_heights  = stage[inlet_indices]-elevation[inlet_indices]
    93         inlet_areas = areas[inlet_indices]
    94 
    95         inlet_water = num.sum(inlet_heights*inlet_areas)
    96 
    97         average_inlet_water = inlet_water/self.exchange_areas[0]
    98 
    99         # Outlet averages
    100         outlet_heights  = stage[outlet_indices]-elevation[outlet_indices]
    101         outlet_areas = areas[outlet_indices]
    102 
    103         outlet_water = num.sum(outlet_heights*outlet_areas)
    104 
    105         average_outlet_water = outlet_water/self.exchange_areas[1]
     169        inlet0 = self.inlets[0]
     170        inlet1 = self.inlets[1]
     171
     172
     173        # Aliases to cell indices
     174        inlet0_indices = inlet0.triangle_indices
     175        inlet1_indices = inlet1.triangle_indices
     176
     177
     178        # Inlet0 averages
     179        inlet0_heights  = self.stage[inlet0_indices]-self.elevation[inlet0_indices]
     180        inlet0_areas    = self.areas[inlet0_indices]
     181
     182        inlet0_water_volume = num.sum(inlet0_heights*inlet0_areas)
     183
     184        average_inlet0_height = inlet0_water_volume/inlet0.area
     185
     186        # Inlet1 averages
     187        inlet1_heights  = self.stage[inlet1_indices]-self.elevation[inlet1_indices]
     188        inlet1_areas    = self.areas[inlet1_indices]
     189
     190        inlet1_water_volume = num.sum(inlet1_heights*inlet1_areas)
     191
     192        average_inlet1_height = inlet1_water_volume/inlet1.area
    106193
    107194
    108195        # Transfer
    109         transfer_water = timestep*inlet_water
    110 
    111         stage[inlet_indices] = elevation[inlet_indices] + average_inlet_water - transfer_water
    112         xmom[inlet_indices] = 0.0
    113         ymom[inlet_indices] = 0.0
    114 
    115 
    116         stage[outlet_indices] = elevation[outlet_indices] + average_outlet_water + transfer_water
    117         xmom[outlet_indices] = 0.0
    118         ymom[outlet_indices] = 0.0
     196        transfer_water = timestep*inlet0_water_volume
     197
     198        self.stage[inlet0_indices] = self.elevation[inlet0_indices] + average_inlet0_height - transfer_water
     199        self.xmom[inlet0_indices] = 0.0
     200        self.ymom[inlet0_indices] = 0.0
     201
     202
     203        self.stage[inlet1_indices] = self.elevation[inlet1_indices] + average_inlet1_height + transfer_water
     204        self.xmom[inlet1_indices] = 0.0
     205        self.ymom[inlet1_indices] = 0.0
    119206
    120207
     
    127214        print self.enquiry_gap_factor
    128215       
    129         for i in [0,1]:
     216        for i, inlet in enumerate(self.inlets):
    130217            print '-------------------------------------'
    131             print 'exchange_region %i' % i
     218            print 'Inlet %i' % i
    132219            print '-------------------------------------'
    133220
    134             print 'exchange triangle indices and centres'
    135             print self.exchange_triangle_indices[i]
    136             print self.domain.get_centroid_coordinates()[self.exchange_triangle_indices[i]]
    137        
    138             print 'end_point'
    139             print self.end_points[i]
    140 
    141 
    142             print 'exchange_polygon'
    143             print self.exchange_polygons[i]
     221            print 'inlet triangle indices and centres'
     222            print inlet.triangle_indices[i]
     223            print self.domain.get_centroid_coordinates()[inlet.triangle_indices[i]]
     224       
     225            print 'polygon'
     226            print inlet.polygon
    144227
    145228            print 'enquiry_point'
    146             print self.enquiry_points[i]
     229            print inlet.enquiry_point
    147230
    148231        print '====================================='
     
    151234 
    152235
    153 
    154     def compute_exchange_triangle_indices(self):
    155 
    156         # Get boundary (in absolute coordinates)
    157         domain = self.domain
    158         bounding_polygon = domain.get_boundary_polygon()
    159         centroids = domain.get_centroid_coordinates(absolute=True)
    160         self.exchange_triangle_indices = []
    161         self.exchange_areas = []
    162 
    163         for i in [0,1]:
    164             exchange_polygon = self.exchange_polygons[i]
    165 
    166             # Check that polygon lies within the mesh.
    167             for point in exchange_polygon:
    168                 msg = 'Point %s in polygon for forcing term' % str(point)
    169                 msg += ' did not fall within the domain boundary.'
    170                 assert is_inside_polygon(point, bounding_polygon), msg
    171 
    172            
    173             exchange_triangle_indices = inside_polygon(centroids, exchange_polygon)
    174 
    175             if len(exchange_triangle_indices) == 0:
    176                 region = 'polygon=%s' % (exchange_polygon)
    177                 msg = 'No triangles have been identified in '
    178                 msg += 'specified region: %s' % region
    179                 raise Exception, msg
    180 
    181             # Compute exchange area as the sum of areas of triangles identified
    182             # by polygon
    183             exchange_area = 0.0
    184             for j in exchange_triangle_indices:
    185                 exchange_area += domain.areas[j]
    186 
    187 
    188             msg = 'Exchange area %f in culvert' % i
    189             msg += ' has area = %f' % exchange_area
    190             assert exchange_area > 0.0
    191 
    192             self.exchange_triangle_indices.append(exchange_triangle_indices)
    193             self.exchange_areas.append(exchange_area)
    194236
    195237
     
    274316
    275317
    276     def compute_enquiry_indices(self):
    277         """Get indices for nearest centroids to self.enquiry_points
    278         """
    279        
    280         domain = self.domain
    281        
    282         enquiry_indices = []                 
    283         for point in self.enquiry_points:
    284             # Find nearest centroid
    285             N = len(domain)   
    286             points = domain.get_centroid_coordinates(absolute=True)
    287 
    288             # Calculate indices in exchange area for this forcing term
    289            
    290             triangle_id = min_dist = sys.maxint
    291             for k in range(N):
    292                 x, y = points[k,:] # Centroid
    293 
    294                 c = point                               
    295                 distance = (x-c[0])**2+(y-c[1])**2
    296                 if distance < min_dist:
    297                     min_dist = distance
    298                     triangle_id = k
    299 
    300                    
    301             if triangle_id < sys.maxint:
    302                 msg = 'found triangle with centroid (%f, %f)'\
    303                     %tuple(points[triangle_id, :])
    304                 msg += ' for point (%f, %f)' %tuple(point)
    305                
    306                 enquiry_indices.append(triangle_id)
    307             else:
    308                 msg = 'Triangle not found for point (%f, %f)' %point
    309                 raise Exception, msg
    310        
    311         self.enquiry_indices = enquiry_indices
    312 
    313        
    314     def check_culvert_inside_domain(self):
    315         """Check that all polygons and enquiry points lie within the mesh.
    316         """
    317         bounding_polygon = self.domain.get_boundary_polygon()
    318         for i in [0, 1]:
    319             for point in list(self.exchange_polygons[i]) + self.enquiry_points:
    320                 msg = 'Point %s did not '\
    321                     %(str(point))
    322                 msg += 'fall within the domain boundary.'
    323                 assert is_inside_polygon(point, bounding_polygon), msg
     318
     319       
     320
    324321           
    325322
  • trunk/anuga_core/source/anuga/structures/testing_culvert.py

    r7962 r7968  
    131131#min_delta_w = sys.maxint
    132132#max_delta_w = -min_delta_w
    133 for t in domain.evolve(yieldstep = 1.0, finaltime = 300):
     133for t in domain.evolve(yieldstep = 1.0, finaltime = 100):
    134134    domain.write_time()
    135135    #delta_w = culvert.inlet.stage - culvert.outlet.stage
Note: See TracChangeset for help on using the changeset viewer.