- Timestamp:
- Feb 27, 2009, 11:54:09 AM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py
r6410 r6428 37 37 triangles = [ [1,0,2], [1,2,3] ] # bac, bce 38 38 39 # Create mesh with two triangles: bac and bce 39 # Create mesh with two triangles: bac and bce 40 40 mesh = Mesh(nodes, triangles) 41 41 … … 56 56 """ 57 57 58 #FIXME: It would be a good idea to use geospatial data as an alternative 59 #input 60 def __init__(self, nodes, triangles, 61 geo_reference=None, 58 # FIXME: It would be a good idea to use geospatial data as an alternative 59 # input 60 def __init__(self, 61 nodes, 62 triangles, 63 geo_reference=None, 62 64 number_of_full_nodes=None, 63 number_of_full_triangles=None, 65 number_of_full_triangles=None, 64 66 verbose=False): 65 67 """Build triangular 2d mesh from nodes and triangle information 66 68 67 69 Input: 68 70 69 71 nodes: x,y coordinates represented as a sequence of 2-tuples or 70 72 a Nx2 numeric array of floats. 71 73 72 74 triangles: sequence of 3-tuples or Mx3 numeric array of 73 75 non-negative integers representing indices into 74 76 the nodes array. 75 77 76 78 georeference (optional): If specified coordinates are 77 79 assumed to be relative to this origin. … … 83 85 In this case it is usefull to specify the number of real (called full) 84 86 nodes and triangles. If omitted they will default to all. 85 86 """ 87 88 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 87 88 """ 89 90 if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain' 89 91 90 92 self.triangles = num.array(triangles, num.int) 91 93 self.nodes = num.array(nodes, num.float) 92 94 93 94 # Register number of elements and nodes 95 # Register number of elements and nodes 95 96 self.number_of_triangles = N = self.triangles.shape[0] 96 self.number_of_nodes = self.nodes.shape[0] 97 98 97 self.number_of_nodes = self.nodes.shape[0] 99 98 100 99 if number_of_full_nodes is None: … … 102 101 else: 103 102 assert int(number_of_full_nodes) 104 self.number_of_full_nodes = number_of_full_nodes 105 103 self.number_of_full_nodes = number_of_full_nodes 106 104 107 105 if number_of_full_triangles is None: 108 self.number_of_full_triangles = self.number_of_triangles 106 self.number_of_full_triangles = self.number_of_triangles 109 107 else: 110 assert int(number_of_full_triangles) 108 assert int(number_of_full_triangles) 111 109 self.number_of_full_triangles = number_of_full_triangles 112 113 114 #print self.number_of_full_nodes, self.number_of_nodes115 #print self.number_of_full_triangles, self.number_of_triangles116 117 118 110 119 111 # FIXME: this stores a geo_reference, but when coords are returned 120 112 # This geo_ref is not taken into account! 121 113 if geo_reference is None: 122 self.geo_reference = Geo_reference() #Use defaults114 self.geo_reference = Geo_reference() # Use defaults 123 115 else: 124 116 self.geo_reference = geo_reference 125 117 126 118 # Input checks 127 msg = 'Triangles must an Mx3 numeric array or a sequence of 3-tuples. '128 msg += 'The supplied array has the shape: %s'\129 % str(self.triangles.shape)119 msg = ('Triangles must an Mx3 numeric array or a sequence of 3-tuples. ' 120 'The supplied array has the shape: %s' 121 % str(self.triangles.shape)) 130 122 assert len(self.triangles.shape) == 2, msg 131 123 132 msg = 'Nodes must an Nx2 numeric array or a sequence of 2-tuples' 133 msg += 'The supplied array has the shape: %s'\ 134 %str(self.nodes.shape) 124 msg = ('Nodes must an Nx2 numeric array or a sequence of 2-tuples' 125 'The supplied array has the shape: %s' % str(self.nodes.shape)) 135 126 assert len(self.nodes.shape) == 2, msg 136 127 … … 138 129 assert max(self.triangles.flat) < self.nodes.shape[0], msg 139 130 140 141 131 # FIXME: Maybe move to statistics? 142 132 # Or use with get_extent 143 xy_extent = [ min(self.nodes[:,0]), min(self.nodes[:,1]),144 max(self.nodes[:,0]), max(self.nodes[:,1])]133 xy_extent = [min(self.nodes[:,0]), min(self.nodes[:,1]), 134 max(self.nodes[:,0]), max(self.nodes[:,1])] 145 135 146 136 self.xy_extent = num.array(xy_extent, num.float) 147 148 137 149 138 # Allocate space for geometric quantities … … 155 144 self.vertex_coordinates = V = self.compute_vertex_coordinates() 156 145 157 158 146 # Initialise each triangle 159 147 if verbose: 160 print 'General_mesh: Computing areas, normals and edgeleng hts'161 148 print 'General_mesh: Computing areas, normals and edgelengths' 149 162 150 for i in range(N): 163 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N)164 151 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' % (i, N) 152 165 153 x0, y0 = V[3*i, :] 166 154 x1, y1 = V[3*i+1, :] 167 x2, y2 = V[3*i+2, :] 155 x2, y2 = V[3*i+2, :] 168 156 169 157 # Area 170 self.areas[i] = abs((x1*y0-x0*y1) +(x2*y1-x1*y2)+(x0*y2-x2*y0))/2171 172 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2)173 msg += ' is degenerate: area == %f' % self.areas[i]158 self.areas[i] = abs((x1*y0-x0*y1) + (x2*y1-x1*y2) + (x0*y2-x2*y0))/2 159 160 msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' % (x0,y0,x1,y1,x2,y2) 161 msg += ' is degenerate: area == %f' % self.areas[i] 174 162 assert self.areas[i] > 0.0, msg 175 176 163 177 164 # Normals … … 184 171 # the first vertex, etc) 185 172 # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle 186 187 n0 = num.array([x2 - x1, y2 - y1], num.float) 173 n0 = num.array([x2-x1, y2-y1], num.float) 188 174 l0 = num.sqrt(num.sum(n0**2)) 189 175 190 n1 = num.array([x0 - x2, y0 -y2], num.float)176 n1 = num.array([x0-x2, y0-y2], num.float) 191 177 l1 = num.sqrt(num.sum(n1**2)) 192 178 193 n2 = num.array([x1 - x0, y1 -y0], num.float)179 n2 = num.array([x1-x0, y1-y0], num.float) 194 180 l2 = num.sqrt(num.sum(n2**2)) 195 181 … … 207 193 self.edgelengths[i, :] = [l0, l1, l2] 208 194 209 210 195 # Build structure listing which trianglse belong to which node. 211 if verbose: print 'Building inverted triangle structure' 196 if verbose: print 'Building inverted triangle structure' 212 197 self.build_inverted_triangle_structure() 213 214 215 198 216 199 def __len__(self): 217 200 return self.number_of_triangles 218 219 201 220 202 def __repr__(self): 221 return 'Mesh: %d vertices, %d triangles'\222 %(self.nodes.shape[0], len(self))203 return ('Mesh: %d vertices, %d triangles' 204 % (self.nodes.shape[0], len(self))) 223 205 224 206 def get_normals(self): 225 207 """Return all normal vectors. 208 226 209 Return normal vectors for all triangles as an Nx6 array 227 210 (ordered as x0, y0, x1, y1, x2, y2 for each triangle) 228 211 """ 212 229 213 return self.normals 230 231 214 232 215 def get_normal(self, i, j): 233 216 """Return normal vector j of the i'th triangle. 217 234 218 Return value is the numeric array slice [x, y] 235 219 """ 220 236 221 return self.normals[i, 2*j:2*j+2] 237 222 238 223 def get_number_of_nodes(self): 239 224 return self.number_of_nodes 240 225 241 226 def get_nodes(self, absolute=False): 242 227 """Return all nodes in mesh. … … 249 234 are to be made absolute by taking georeference into account 250 235 Default is False as many parts of ANUGA expects relative coordinates. 251 (To see which, switch to default absolute=True and run tests). 236 (To see which, switch to default absolute=True and run tests). 252 237 """ 253 238 … … 257 242 if not self.geo_reference.is_absolute(): 258 243 V = self.geo_reference.get_absolute(V) 259 244 260 245 return V 261 262 def get_node(self, i, 263 absolute=False): 246 247 def get_node(self, i, absolute=False): 264 248 """Return node coordinates for triangle i. 265 249 … … 267 251 are to be made absolute by taking georeference into account 268 252 Default is False as many parts of ANUGA expects relative coordinates. 269 (To see which, switch to default absolute=True and run tests). 270 """ 271 272 253 (To see which, switch to default absolute=True and run tests). 254 """ 255 273 256 V = self.nodes[i,:] 274 257 if absolute is True: 275 258 if not self.geo_reference.is_absolute(): 276 return V + num.array([self.geo_reference.get_xllcorner(), 277 self.geo_reference.get_yllcorner()], num.float) 278 else: 279 return V 280 else: 281 return V 282 283 284 285 def get_vertex_coordinates(self, 286 triangle_id=None, 287 absolute=False): 288 """Return vertex coordinates for all triangles. 289 259 V += num.array([self.geo_reference.get_xllcorner(), 260 self.geo_reference.get_yllcorner()], num.float) 261 return V 262 263 def get_vertex_coordinates(self, triangle_id=None, absolute=False): 264 """Return vertex coordinates for all triangles. 265 290 266 Return all vertex coordinates for all triangles as a 3*M x 2 array 291 267 where the jth vertex of the ith triangle is located in row 3*i+j and … … 294 270 if triangle_id is specified (an integer) the 3 vertex coordinates 295 271 for triangle_id are returned. 296 272 297 273 Boolean keyword argument absolute determines whether coordinates 298 274 are to be made absolute by taking georeference into account 299 275 Default is False as many parts of ANUGA expects relative coordinates. 300 276 """ 301 277 302 278 V = self.vertex_coordinates 303 279 304 if triangle_id is None: 280 if triangle_id is None: 305 281 if absolute is True: 306 282 if not self.geo_reference.is_absolute(): 307 283 V = self.geo_reference.get_absolute(V) 308 309 284 return V 310 285 else: … … 313 288 assert int(i) == i, msg 314 289 assert 0 <= i < self.number_of_triangles 315 316 i3 = 3*i 290 291 i3 = 3*i 317 292 if absolute is True and not self.geo_reference.is_absolute(): 318 293 offset=num.array([self.geo_reference.get_xllcorner(), … … 323 298 else: 324 299 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float) 325 326 327 300 328 301 def get_vertex_coordinate(self, i, j, absolute=False): … … 333 306 msg = 'vertex id j must be an integer in [0,1,2]' 334 307 assert j in [0,1,2], msg 335 336 V = self.get_vertex_coordinates(triangle_id=i, 337 absolute=absolute) 308 309 V = self.get_vertex_coordinates(triangle_id=i, absolute=absolute) 338 310 return V[j,:] 339 340 341 311 342 312 def compute_vertex_coordinates(self): … … 358 328 return vertex_coordinates 359 329 360 361 362 330 def get_triangles(self, indices=None): 363 331 """Get mesh triangles. … … 375 343 if indices is None: 376 344 return self.triangles 377 #indices = range(M)378 345 379 346 return num.take(self.triangles, indices, axis=0) 380 381 382 347 383 348 def get_disconnected_triangles(self): … … 398 363 399 364 The triangles created will have the format 400 401 365 [[0,1,2], 402 366 [3,4,5], 403 367 [6,7,8], 404 368 ... 405 [3*M-3 3*M-2 3*M-1]] 369 [3*M-3 3*M-2 3*M-1]] 406 370 """ 407 371 408 372 M = len(self) # Number of triangles 409 373 K = 3*M # Total number of unique vertices 410 411 T = num.reshape(num.arange(K, dtype=num.int), (M,3)) 412 413 return T 414 415 416 417 def get_unique_vertices(self, indices=None): 374 return num.reshape(num.arange(K, dtype=num.int), (M,3)) 375 376 def get_unique_vertices(self, indices=None): 418 377 """FIXME(Ole): This function needs a docstring""" 419 378 … … 421 380 unique_verts = {} 422 381 for triangle in triangles: 423 #print 'triangle(%s)=%s' % (type(triangle), str(triangle))424 382 unique_verts[triangle[0]] = 0 425 383 unique_verts[triangle[1]] = 0 … … 427 385 return unique_verts.keys() 428 386 429 430 387 def get_triangles_and_vertices_per_node(self, node=None): 431 388 """Get triangles associated with given node. … … 441 398 # Get index for this node 442 399 first = num.sum(self.number_of_triangles_per_node[:node]) 443 400 444 401 # Get number of triangles for this node 445 402 count = self.number_of_triangles_per_node[node] … … 459 416 # working directly with the inverted triangle structure 460 417 for i in range(self.number_of_full_nodes): 461 462 418 L = self.get_triangles_and_vertices_per_node(node=i) 463 464 419 triangle_list.append(L) 465 420 466 421 return triangle_list 467 468 422 469 423 def build_inverted_triangle_structure(self): … … 475 429 listing for each node how many triangles use it. N is the number of 476 430 nodes in mesh. 477 431 478 432 vertex_value_indices: An array of length M listing indices into 479 433 triangles ordered by node number. The (triangle_id, vertex_id) … … 483 437 used to average vertex values efficiently. 484 438 485 486 439 Example: 487 488 440 a = [0.0, 0.0] # node 0 489 441 b = [0.0, 2.0] # node 1 … … 492 444 e = [2.0, 2.0] # node 4 493 445 f = [4.0, 0.0] # node 5 494 495 446 nodes = array([a, b, c, d, e, f]) 496 497 #bac, bce, ecf, dbe, daf, dae 498 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 499 500 501 For this structure 502 447 448 # bac, bce, ecf, dbe 449 triangles = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]]) 450 451 For this structure: 503 452 number_of_triangles_per_node = [1 3 3 1 3 1] 504 453 which means that node a has 1 triangle associated with it, node b 505 454 has 3, node has 3 and so on. 506 455 507 456 vertex_value_indices = [ 1 0 3 10 2 4 7 9 5 6 11 8] 508 457 which reflects the fact that … … 518 467 and by triangle 2, vertex 0 (index = 6) 519 468 and by triangle 3, vertex 2 (index = 11) 520 node 5 is used by triangle 2, vertex 2 (index = 8) 521 469 node 5 is used by triangle 2, vertex 2 (index = 8) 522 470 523 471 Preconditions: … … 526 474 Postcondition: 527 475 self.number_of_triangles_per_node is built 528 self.vertex_value_indices is built 476 self.vertex_value_indices is built 529 477 """ 530 478 … … 541 489 542 490 # Register (triangle, vertex) indices for each node 543 vertexlist = [None] *self.number_of_full_nodes491 vertexlist = [None] * self.number_of_full_nodes 544 492 for volume_id in range(self.number_of_full_triangles): 545 546 493 a = self.triangles[volume_id, 0] 547 494 b = self.triangles[volume_id, 1] 548 495 c = self.triangles[volume_id, 2] 549 496 550 for vertex_id, node_id in enumerate([a, b,c]):497 for vertex_id, node_id in enumerate([a, b, c]): 551 498 if vertexlist[node_id] is None: 552 499 vertexlist[node_id] = [] 553 554 vertexlist[node_id].append( (volume_id, vertex_id) ) 555 500 vertexlist[node_id].append((volume_id, vertex_id)) 556 501 557 502 # Build inverted triangle index array … … 561 506 for volume_id, vertex_id in vertices: 562 507 vertex_value_indices[k] = 3*volume_id + vertex_id 563 564 508 k += 1 565 509 … … 568 512 self.vertex_value_indices = vertex_value_indices 569 513 570 571 572 573 514 def get_extent(self, absolute=False): 574 515 """Return min and max of all x and y coordinates … … 577 518 are to be made absolute by taking georeference into account 578 519 """ 579 580 581 520 582 521 C = self.get_vertex_coordinates(absolute=absolute) … … 588 527 ymin = min(Y.flat) 589 528 ymax = max(Y.flat) 590 #print "C",C 529 591 530 return xmin, xmax, ymin, ymax 592 531 593 532 def get_areas(self): 594 """Get areas of all individual triangles.595 """ 596 return self.areas 533 """Get areas of all individual triangles.""" 534 535 return self.areas 597 536 598 537 def get_area(self): 599 """Return total area of mesh 600 """ 538 """Return total area of mesh""" 601 539 602 540 return num.sum(self.areas) … … 604 542 def set_georeference(self, g): 605 543 self.geo_reference = g 606 544 607 545 def get_georeference(self): 608 546 return self.geo_reference 609 610 611 547
Note: See TracChangeset
for help on using the changeset viewer.
