Changeset 2705
- Timestamp:
- Apr 12, 2006, 10:58:17 AM (19 years ago)
- Location:
- development/pyvolution-1d
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
development/pyvolution-1d/domain.py
r2229 r2705 7 7 Geoscience Australia 8 8 """ 9 from generic_boundary_conditions import * 9 10 10 11 class Domain: 11 12 12 def __init__(self, coordinates, conserved_quantities = None, 13 other_quantities = None): 13 def __init__(self, coordinates, boundary = None, 14 conserved_quantities = None, 15 other_quantities = None, 16 tagged_elements = None): 14 17 """ 15 18 Build 1D elements from x coordinates … … 25 28 self.number_of_elements = N = len(self.coordinates)-1 26 29 27 30 #Allocate space for neighbour and boundary structures 31 self.neighbours = zeros((N, 2), Int) 32 self.neighbour_edges = zeros((N, 2), Int) 33 self.number_of_boundaries = zeros(N, Int) 34 self.surrogate_neighbours = zeros((N, 2), Int) 35 28 36 #Allocate space for geometric quantities 29 37 self.vertices = zeros((N, 2), Float) … … 43 51 44 52 self.areas[i] = (xr-xl) 45 53 46 54 ## print 'N', N 47 55 ## print 'Centroid', self.centroids 48 56 ## print 'Areas', self.areas 49 57 ## print 'Vertex_Coordinates', self.vertices 58 59 #Initialise Neighbours (-1 means that it is a boundary neighbour) 60 self.neighbours[i, :] = [-1, -1] 61 #Initialise edge ids of neighbours 62 #In case of boundaries this slot is not used 63 self.neighbour_edges[i, :] = [-1, -1] 64 65 #Build neighbour structure 66 self.build_neighbour_structure() 67 68 #Build surrogate neighbour structure 69 self.build_surrogate_neighbour_structure() 70 71 #Build boundary dictionary mapping (id, edge) to symbolic tags 72 self.build_boundary_dictionary(boundary) 73 74 #Build tagged element dictionary mapping (tag) to array of elements 75 self.build_tagged_elements_dictionary(tagged_elements) 76 77 78 79 from quantity import Quantity, Conserved_quantity 80 81 #List of quantity names entering 82 #the conservation equations 83 #(Must be a subset of quantities) 84 if conserved_quantities is None: 85 self.conserved_quantities = [] 86 else: 87 self.conserved_quantities = conserved_quantities 88 89 if other_quantities is None: 90 self.other_quantities = [] 91 else: 92 self.other_quantities = other_quantities 93 94 95 #Build dictionary of Quantity instances keyed by quantity names 96 self.quantities = {} 97 98 #FIXME: remove later - maybe OK, though.... 99 for name in self.conserved_quantities: 100 self.quantities[name] = Conserved_quantity(self) 101 for name in self.other_quantities: 102 self.quantities[name] = Quantity(self) 103 104 #Create an empty list for explicit forcing terms 105 self.forcing_terms = [] 106 107 def build_neighbour_structure(self): 108 """Update all registered triangles to point to their neighbours. 109 110 Also, keep a tally of the number of boundaries for each triangle 111 112 Postconditions: 113 neighbours and neighbour_edges is populated 114 number_of_boundaries integer array is defined. 115 """ 116 117 #Step 1: 118 #Build dictionary mapping from segments (2-tuple of points) 119 #to left hand side edge (facing neighbouring triangle) 120 121 N = self.number_of_elements 122 neighbourdict = {} 123 l_edge = 0 124 r_edge = 1 125 for i in range(N): 126 127 #Register all segments as keys mapping to current triangle 128 #and segment id 129 #a = self.triangles[i, 0] 130 #b = self.triangles[i, 1] 131 #c = self.triangles[i, 2] 132 a = self.vertices[i,0] 133 b = self.vertices[i,1] 134 135 """ 136 if neighbourdict.has_key((a,b)): 137 msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0]) 138 raise msg 139 if neighbourdict.has_key((b,c)): 140 msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0]) 141 raise msg 142 if neighbourdict.has_key((c,a)): 143 msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0]) 144 raise msg 145 """ 146 #neighbourdict[a,b] = (i, 2) #(id, edge) 147 #neighbourdict[b,c] = (i, 0) #(id, edge) 148 #neighbourdict[c,a] = (i, 1) #(id, edge) 149 #neighbourdict[a,b] = (i, 1) #(id, edge) 150 #neighbourdict[b,a] = (i, 0) #(id, edge) 151 neighbourdict[a,l_edge] = (i, 0) #(id, edge) 152 neighbourdict[b,r_edge] = (i, 1) #(id, edge) 153 154 155 #Step 2: 156 #Go through triangles again, but this time 157 #reverse direction of segments and lookup neighbours. 158 for i in range(N): 159 #a = self.triangles[i, 0] 160 #b = self.triangles[i, 1] 161 #c = self.triangles[i, 2] 162 163 a = self.vertices[i,0] 164 b = self.vertices[i,1] 165 166 #self.number_of_boundaries[i] = 3 167 self.number_of_boundaries[i] = 2 168 if neighbourdict.has_key((b,l_edge)): 169 self.neighbours[i, 1] = neighbourdict[b,l_edge][0] 170 self.neighbour_edges[i, 1] = neighbourdict[b,l_edge][1] 171 self.number_of_boundaries[i] -= 1 172 173 if neighbourdict.has_key((a,r_edge)): 174 self.neighbours[i, 0] = neighbourdict[a,r_edge][0] 175 self.neighbour_edges[i, 0] = neighbourdict[a,r_edge][1] 176 self.number_of_boundaries[i] -= 1 177 #if neighbourdict.has_key((b,a)): 178 # self.neighbours[i, 1] = neighbourdict[b,a][0] 179 # self.neighbour_edges[i, 1] = neighbourdict[b,a][1] 180 # self.number_of_boundaries[i] -= 1 181 182 #if neighbourdict.has_key((c,b)): 183 # self.neighbours[i, 0] = neighbourdict[c,b][0] 184 # self.neighbour_edges[i, 0] = neighbourdict[c,b][1] 185 # self.number_of_boundaries[i] -= 1 186 187 #if neighbourdict.has_key((a,b)): 188 # self.neighbours[i, 0] = neighbourdict[a,b][0] 189 # self.neighbour_edges[i, 0] = neighbourdict[a,b][1] 190 # self.number_of_boundaries[i] -= 1 191 192 def build_surrogate_neighbour_structure(self): 193 """Build structure where each triangle edge points to its neighbours 194 if they exist. Otherwise point to the triangle itself. 195 196 The surrogate neighbour structure is useful for computing gradients 197 based on centroid values of neighbours. 198 199 Precondition: Neighbour structure is defined 200 Postcondition: 201 Surrogate neighbour structure is defined: 202 surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to 203 triangles. 204 205 """ 206 207 N = self.number_of_elements 208 for i in range(N): 209 #Find all neighbouring volumes that are not boundaries 210 #for k in range(3): 211 for k in range(2): 212 if self.neighbours[i, k] < 0: 213 self.surrogate_neighbours[i, k] = i #Point this triangle 214 else: 215 self.surrogate_neighbours[i, k] = self.neighbours[i, k] 216 217 def build_boundary_dictionary(self, boundary = None): 218 """Build or check the dictionary of boundary tags. 219 self.boundary is a dictionary of tags, 220 keyed by volume id and edge: 221 { (id, edge): tag, ... } 222 223 Postconditions: 224 self.boundary is defined. 225 """ 226 227 from config import default_boundary_tag 228 229 if boundary is None: 230 boundary = {} 231 for vol_id in range(self.number_of_elements): 232 #for edge_id in range(0, 3): 233 for edge_id in range(0, 2): 234 if self.neighbours[vol_id, edge_id] < 0: 235 boundary[(vol_id, edge_id)] = default_boundary_tag 236 else: 237 #Check that all keys in given boundary exist 238 for vol_id, edge_id in boundary.keys(): 239 msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id) 240 a, b = self.neighbours.shape 241 assert vol_id < a and edge_id < b, msg 242 243 #FIXME: This assert violates internal boundaries (delete it) 244 #msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id) 245 #assert self.neighbours[vol_id, edge_id] < 0, msg 246 247 #Check that all boundary segments are assigned a tag 248 for vol_id in range(self.number_of_elements): 249 #for edge_id in range(0, 3): 250 for edge_id in range(0, 2): 251 if self.neighbours[vol_id, edge_id] < 0: 252 if not boundary.has_key( (vol_id, edge_id) ): 253 msg = 'WARNING: Given boundary does not contain ' 254 msg += 'tags for edge (%d, %d). '\ 255 %(vol_id, edge_id) 256 msg += 'Assigning default tag (%s).'\ 257 %default_boundary_tag 258 259 #FIXME: Print only as per verbosity 260 #print msg 261 262 #FIXME: Make this situation an error in the future 263 #and make another function which will 264 #enable default boundary-tags where 265 #tags a not specified 266 boundary[ (vol_id, edge_id) ] =\ 267 default_boundary_tag 268 269 270 271 self.boundary = boundary 272 273 def build_tagged_elements_dictionary(self, tagged_elements = None): 274 """Build the dictionary of element tags. 275 self.tagged_elements is a dictionary of element arrays, 276 keyed by tag: 277 { (tag): [e1, e2, e3..] } 278 279 Postconditions: 280 self.element_tag is defined 281 """ 282 from Numeric import array, Int 283 284 if tagged_elements is None: 285 tagged_elements = {} 286 else: 287 #Check that all keys in given boundary exist 288 for tag in tagged_elements.keys(): 289 tagged_elements[tag] = array(tagged_elements[tag]).astype(Int) 290 291 msg = 'Not all elements exist. ' 292 assert max(tagged_elements[tag]) < self.number_of_elements, msg 293 #print "tagged_elements", tagged_elements 294 self.tagged_elements = tagged_elements 295 296 def get_boundary_tags(self): 297 """Return list of available boundary tags 298 """ 299 300 tags = {} 301 for v in self.boundary.values(): 302 tags[v] = 1 303 304 return tags.keys() 305 306 307 def get_conserved_quantities(self, vol_id, vertex=None, edge=None): 308 """Get conserved quantities at volume vol_id 309 310 If vertex is specified use it as index for vertex values 311 If edge is specified use it as index for edge values 312 If neither are specified use centroid values 313 If both are specified an exeception is raised 314 315 Return value: Vector of length == number_of_conserved quantities 316 317 """ 318 319 from Numeric import zeros, Float 320 321 if not (vertex is None or edge is None): 322 msg = 'Values for both vertex and edge was specified.' 323 msg += 'Only one (or none) is allowed.' 324 raise msg 325 326 q = zeros( len(self.conserved_quantities), Float) 327 328 for i, name in enumerate(self.conserved_quantities): 329 Q = self.quantities[name] 330 if vertex is not None: 331 q[i] = Q.vertex_values[vol_id, vertex] 332 elif edge is not None: 333 q[i] = Q.edge_values[vol_id, edge] 334 else: 335 q[i] = Q.centroid_values[vol_id] 336 337 return q 338 50 339 51 340 def get_centroids(self): … … 79 368 80 369 return self.areas[elem_id] 370 371 def set_quantity(self, name, *args, **kwargs): 372 """Set values for named quantity 373 374 375 One keyword argument is documented here: 376 expression = None, # Arbitrary expression 377 378 expression: 379 Arbitrary expression involving quantity names 380 381 See Quantity.set_values for further documentation. 382 """ 383 384 #FIXME (Ole): Allow new quantities here 385 #from quantity import Quantity, Conserved_quantity 386 #Create appropriate quantity object 387 ##if name in self.conserved_quantities: 388 ## self.quantities[name] = Conserved_quantity(self) 389 ##else: 390 ## self.quantities[name] = Quantity(self) 391 392 393 #Do the expression stuff 394 if kwargs.has_key('expression'): 395 expression = kwargs['expression'] 396 del kwargs['expression'] 397 398 Q = self.create_quantity_from_expression(expression) 399 kwargs['quantity'] = Q 400 401 402 #Assign values 403 self.quantities[name].set_values(*args, **kwargs) 404 405 def set_boundary(self, boundary_map): 406 """Associate boundary objects with tagged boundary segments. 407 408 Input boundary_map is a dictionary of boundary objects keyed 409 by symbolic tags to matched against tags in the internal dictionary 410 self.boundary. 411 412 As result one pointer to a boundary object is stored for each vertex 413 in the list self.boundary_objects. 414 More entries may point to the same boundary object 415 416 Schematically the mapping is from two dictionaries to one list 417 where the index is used as pointer to the boundary_values arrays 418 within each quantity. 419 420 self.boundary: (vol_id, edge_id): tag 421 boundary_map (input): tag: boundary_object 422 ---------------------------------------------- 423 self.boundary_objects: ((vol_id, edge_id), boundary_object) 424 425 426 Pre-condition: 427 self.boundary has been built. 428 429 Post-condition: 430 self.boundary_objects is built 431 432 If a tag from the domain doesn't appear in the input dictionary an 433 exception is raised. 434 However, if a tag is not used to the domain, no error is thrown. 435 FIXME: This would lead to implementation of a 436 default boundary condition 437 438 Note: If a segment is listed in the boundary dictionary and if it is 439 not None, it *will* become a boundary - 440 even if there is a neighbouring triangle. 441 This would be the case for internal boundaries 442 443 Boundary objects that are None will be skipped. 444 445 FIXME: If set_boundary is called multiple times and if Boundary 446 object is changed into None, the neighbour structure will not be 447 restored!!! 448 """ 449 450 self.boundary_objects = [] 451 452 453 454 455 456 self.boundary_map = boundary_map #Store for use with eg. boundary_stats. 457 458 #FIXME: Try to remove the sorting and fix test_mesh.py 459 x = self.boundary.keys() 460 x.sort() 461 462 #Loop through edges that lie on the boundary and associate them with 463 #callable boundary objects depending on their tags 464 for k, (vol_id, edge_id) in enumerate(x): 465 tag = self.boundary[ (vol_id, edge_id) ] 466 467 if boundary_map.has_key(tag): 468 B = boundary_map[tag] #Get callable boundary object 469 470 if B is not None: 471 self.boundary_objects.append( ((vol_id, edge_id), B) ) 472 self.neighbours[vol_id, edge_id] = -len(self.boundary_objects) 473 else: 474 pass 475 #FIXME: Check and perhaps fix neighbour structure 476 477 else: 478 msg = 'ERROR (domain.py): Tag "%s" has not been ' %tag 479 msg += 'bound to a boundary object.\n' 480 msg += 'All boundary tags defined in domain must appear ' 481 msg += 'in the supplied dictionary.\n' 482 msg += 'The tags are: %s' %self.get_boundary_tags() 483 raise msg 484 485 486 487 def check_integrity(self): 488 #Mesh.check_integrity(self) 489 490 for quantity in self.conserved_quantities: 491 msg = 'Conserved quantities must be a subset of all quantities' 492 assert quantity in self.quantities, msg 493 494 ##assert hasattr(self, 'boundary_objects') 495 496 def update_boundary(self): 497 """Go through list of boundary objects and update boundary values 498 for all conserved quantities on boundary. 499 """ 500 501 #FIXME: Update only those that change (if that can be worked out) 502 #FIXME: Boundary objects should not include ghost nodes. 503 for i, ((vol_id, edge_id), B) in enumerate(self.boundary_objects): 504 q = B.evaluate(vol_id, edge_id) 505 506 for j, name in enumerate(self.conserved_quantities): 507 Q = self.quantities[name] 508 Q.boundary_values[i] = q[j] 81 509 82 510 … … 101 529 #D2 = Domain(points2) 102 530 103 104 105 106 107
Note: See TracChangeset
for help on using the changeset viewer.