Changeset 2709 for inundation/pyvolution/mesh.py
 Timestamp:
 Apr 12, 2006, 2:28:05 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/pyvolution/mesh.py
r2683 r2709 6 6 7 7 from general_mesh import General_mesh 8 from math import pi, sqrt 9 8 10 9 11 class Mesh(General_mesh): … … 125 127 #of inscribed circle is computed 126 128 127 from math import sqrt128 129 a = sqrt((x0x1)**2+(y0y1)**2) 129 130 b = sqrt((x1x2)**2+(y1y2)**2) … … 165 166 def set_to_inscribed_circle(self,safety_factor = 1): 166 167 #FIXME phase out eventually 167 from math import sqrt168 168 N = self.number_of_elements 169 169 V = self.vertex_coordinates … … 414 414 415 415 416 def get_boundary_polygon(self): 417 """Return bounding polygon as a list of points 418 419 FIXME: If triangles are listed as discontinuous 420 (e.g vertex coordinates listed multiple times), 421 this may not work as expected. 422 """ 416 def get_boundary_polygon(self, verbose = False): 417 """Return bounding polygon for mesh (counter clockwise) 418 419 Using the mesh boundary, derive a bounding polygon for this mesh. 420 If multiple vertex values are present, the algorithm will select the 421 path that contains the mesh. 422 """ 423 423 424 from Numeric import allclose, sqrt, array, minimum, maximum 424 425 426 427 #V = self.get_vertex_coordinates() 425 from utilities.numerical_tools import angle, ensure_numeric 426 427 428 # Get mesh extent 429 xmin, xmax, ymin, ymax = self.get_extent() 430 pmin = ensure_numeric([xmin, ymin]) 431 pmax = ensure_numeric([xmax, ymax]) 432 433 434 # Assemble dictionary of boundary segments and choose starting point 428 435 segments = {} 429 430 #pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1])) 431 #pmax = (max(self.coordinates[:,0]), max(self.coordinates[:,1])) 432 433 #FIXME:Can this be written more compactly, e.g. 434 #using minimum and maximium? 435 pmin = array( [min(self.coordinates[:,0]), 436 min(self.coordinates[:,1]) ] ) 437 438 pmax = array( [max(self.coordinates[:,0]), 439 max(self.coordinates[:,1]) ] ) 440 441 mindist = sqrt(sum( (pmaxpmin)**2 )) 436 inverse_segments = {} 437 mindist = sqrt(sum((pmaxpmin)**2)) #Start value across entire mesh 442 438 for i, edge_id in self.boundary.keys(): 443 # Find vertex ids for boundary segment439 # Find vertex ids for boundary segment 444 440 if edge_id == 0: a = 1; b = 2 445 441 if edge_id == 1: a = 2; b = 0 446 442 if edge_id == 2: a = 0; b = 1 447 443 448 A = tuple(self.get_vertex_coordinate(i, a))449 B = tuple(self.get_vertex_coordinate(i, b))450 451 # Take the point closest to pmin as starting point452 # Note: Could be arbitrary, but nice to have453 # a unique way of selecting454 dist_A = sqrt(sum( (Apmin)**2))455 dist_B = sqrt(sum( (Bpmin)**2))456 457 #Find minimalpoint444 A = self.get_vertex_coordinate(i, a) # Start 445 B = self.get_vertex_coordinate(i, b) # End 446 447 # Take the point closest to pmin as starting point 448 # Note: Could be arbitrary, but nice to have 449 # a unique way of selecting 450 dist_A = sqrt(sum((Apmin)**2)) 451 dist_B = sqrt(sum((Bpmin)**2)) 452 453 #Find lower leftmost point 458 454 if dist_A < mindist: 459 455 mindist = dist_A … … 464 460 465 461 462 # Sanity check 466 463 if p0 is None: 467 raise 'Weird' 468 p0 = A #We need a starting point (FIXME) 469 470 print 'A', A 471 print 'B', B 472 print 'pmin', pmin 473 print 474 475 segments[A] = B 464 raise Exception('Impossible') 465 466 467 # Register potential paths from A to B 468 if not segments.has_key(tuple(A)): 469 segments[tuple(A)] = [] # Empty list for candidate points 470 471 segments[tuple(A)].append(B) 472 473 476 474 477 475 478 476 #Start with smallest point and follow boundary (counter clock wise) 479 polygon = [p0] 480 while len(polygon) < len(self.boundary): 481 p1 = segments[p0] 477 polygon = [p0] # Storage for final boundary polygon 478 point_registry = {} # Keep track of storage to avoid multiple runs around boundary 479 # This will only be the case if there are more than one candidate 480 # FIXME (Ole): Perhaps we can do away with polygon and use 481 # only point_registry to save space. 482 483 point_registry[tuple(p0)] = 0 484 485 #while len(polygon) < len(self.boundary): 486 while len(point_registry) < len(self.boundary): 487 488 candidate_list = segments[tuple(p0)] 489 if len(candidate_list) > 1: 490 # Multiple points detected 491 # Take the candidate that is furthest to the clockwise direction, 492 # as that will follow the boundary. 493 494 495 if verbose: 496 print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list) 497 498 499 # Choose vector against which all angles will be measured 500 if len(polygon) > 1: 501 v_prev = p0  polygon[2] # Vector that leads to p0 502 else: 503 # FIXME (Ole): What do we do if the first point has multiple candidates? 504 # Being the lower left corner, perhaps we can use the 505 # vector [1, 0], but I really don't know if this is completely 506 # watertight. 507 # Another option might be v_prev = [1.0, 0.0] 508 v_prev = [1.0, 0.0] 509 510 511 512 # Choose candidate with minimum angle 513 minimum_angle = 2*pi 514 for pc in candidate_list: 515 vc = pcp0 # Candidate vector 516 517 # Angle between each candidate and the previous vector in [pi, pi] 518 ac = angle(vc, v_prev) 519 if ac > pi: ac = piac 520 521 # take the minimal angle corresponding to the rightmost vector 522 if ac < minimum_angle: 523 minimum_angle = ac 524 p1 = pc # Best candidate 525 526 527 if verbose is True: 528 print ' Best candidate %s, angle %f' %(p1, minimum_angle*180/pi) 529 530 else: 531 p1 = candidate_list[0] 532 533 if point_registry.has_key(tuple(p1)): 534 # We have completed the boundary polygon  yeehaa 535 break 536 else: 537 point_registry[tuple(p1)] = len(point_registry) 538 482 539 polygon.append(p1) 483 540 p0 = p1 541 484 542 485 543 return polygon … … 494 552 495 553 from config import epsilon 496 from math import pi497 554 from utilities.numerical_tools import anglediff 498 555 … … 537 594 538 595 #Normalise 539 from math import sqrt540 596 l_u = sqrt(u[0]*u[0] + u[1]*u[1]) 541 597 l_v = sqrt(v[0]*v[0] + v[1]*v[1])
Note: See TracChangeset
for help on using the changeset viewer.