Changeset 7968
 Timestamp:
 Aug 23, 2010, 2:48:30 PM (9 years ago)
 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 884 884 return i 885 885 886 return 886 msg = 'Point %s not found within a triangle' %str(point) 887 raise Exception, msg 888 887 889 888 890 
trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7962 r7968 20 20 class Above_interval(Exception): pass 21 21 22 22 23 24 class 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 23 92 24 93 class Generic_box_culvert: … … 60 129 self.filename = None 61 130 62 # Create the fundamental culvert polygons from polygon131 # Create the fundamental culvert polygons and create inlet objects 63 132 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() 70 158 71 159 … … 79 167 80 168 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 106 193 107 194 108 195 # Transfer 109 transfer_water = timestep*inlet _water110 111 s tage[inlet_indices] = elevation[inlet_indices] + average_inlet_water transfer_water112 xmom[inlet_indices]= 0.0113 ymom[inlet_indices]= 0.0114 115 116 s tage[outlet_indices] = elevation[outlet_indices] + average_outlet_water+ transfer_water117 xmom[outlet_indices]= 0.0118 ymom[outlet_indices]= 0.0196 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 119 206 120 207 … … 127 214 print self.enquiry_gap_factor 128 215 129 for i in [0,1]:216 for i, inlet in enumerate(self.inlets): 130 217 print '' 131 print ' exchange_region%i' % i218 print 'Inlet %i' % i 132 219 print '' 133 220 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 144 227 145 228 print 'enquiry_point' 146 print self.enquiry_points[i]229 print inlet.enquiry_point 147 230 148 231 print '=====================================' … … 151 234 152 235 153 154 def compute_exchange_triangle_indices(self):155 156 # Get boundary (in absolute coordinates)157 domain = self.domain158 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), msg171 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' % region179 raise Exception, msg180 181 # Compute exchange area as the sum of areas of triangles identified182 # by polygon183 exchange_area = 0.0184 for j in exchange_triangle_indices:185 exchange_area += domain.areas[j]186 187 188 msg = 'Exchange area %f in culvert' % i189 msg += ' has area = %f' % exchange_area190 assert exchange_area > 0.0191 192 self.exchange_triangle_indices.append(exchange_triangle_indices)193 self.exchange_areas.append(exchange_area)194 236 195 237 … … 274 316 275 317 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 = (xc[0])**2+(yc[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 324 321 325 322 
trunk/anuga_core/source/anuga/structures/testing_culvert.py
r7962 r7968 131 131 #min_delta_w = sys.maxint 132 132 #max_delta_w = min_delta_w 133 for t in domain.evolve(yieldstep = 1.0, finaltime = 300):133 for t in domain.evolve(yieldstep = 1.0, finaltime = 100): 134 134 domain.write_time() 135 135 #delta_w = culvert.inlet.stage  culvert.outlet.stage
Note: See TracChangeset
for help on using the changeset viewer.