Changeset 4515
 Timestamp:
 May 30, 2007, 5:25:51 AM (17 years ago)
 Location:
 anuga_core/source/anuga/abstract_2d_finite_volumes
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r4475 r4515 444 444 445 445 Using the mesh boundary, derive a bounding polygon for this mesh. 446 If multiple vertex values are present , the algorithm will select the447 path that contains the mesh.446 If multiple vertex values are present (vertices stored uniquely), 447 the algorithm will select the path that contains the entire mesh. 448 448 449 449 All points are in absolute UTM coordinates … … 451 451 452 452 from Numeric import allclose, sqrt, array, minimum, maximum 453 from anuga.utilities.numerical_tools import angle, ensure_numeric 453 from anuga.utilities.numerical_tools import angle, ensure_numeric 454 454 455 455 … … 464 464 inverse_segments = {} 465 465 p0 = None 466 mindist = sqrt(sum((pmaxpmin)**2)) # Start value across entire mesh466 mindist = sqrt(sum((pmaxpmin)**2)) # Start value across entire mesh 467 467 for i, edge_id in self.boundary.keys(): 468 468 # Find vertex ids for boundary segment … … 481 481 dist_B = sqrt(sum((Bpmin)**2)) 482 482 483 # Find lower leftmost point483 # Find lower leftmost point 484 484 if dist_A < mindist: 485 485 mindist = dist_A … … 497 497 # Register potential paths from A to B 498 498 if not segments.has_key(tuple(A)): 499 segments[tuple(A)] = [] # Empty list for candidate points 500 499 segments[tuple(A)] = [] # Empty list for candidate points 500 501 501 segments[tuple(A)].append(B) 502 502 503 503 504 # Start with smallest point and follow boundary (counter clock wise)504 # Start with smallest point and follow boundary (counter clock wise) 505 505 polygon = [p0] # Storage for final boundary polygon 506 point_registry = {} # Keep track of storage to avoid multiple runs around507 # boundary. This will only be the case if there are 508 #more than one candidate.506 point_registry = {} # Keep track of storage to avoid multiple runs 507 # around boundary. This will only be the case if 508 # there are more than one candidate. 509 509 # FIXME (Ole): Perhaps we can do away with polygon 510 510 # and use only point_registry to save space. … … 512 512 point_registry[tuple(p0)] = 0 513 513 514 #while len(polygon) < len(self.boundary):515 514 while len(point_registry) < len(self.boundary): 516 515 517 516 candidate_list = segments[tuple(p0)] 518 517 if len(candidate_list) > 1: 519 # Multiple points detected (this will be the case for meshes with 520 # duplicate points as those used for discontinuous triangles). 521 # Take the candidate that is furthest to the clockwise direction, 522 # as that will follow the boundary. 523 518 # Multiple points detected (this will be the case for meshes 519 # with duplicate points as those used for discontinuous 520 # triangles with vertices stored uniquely). 521 # Take the candidate that is furthest to the clockwise 522 # direction, as that will follow the boundary. 523 # 524 # This will also be the case for pathological triangles 525 # that have no neighbours. 524 526 525 527 if verbose: … … 527 529 %(str(p0), candidate_list) 528 530 529 530 531 # Check that previous are not in candidate list 531 532 #for p in candidate_list: 532 533 # assert not allclose(p0, p) 533 534 534 535 535 # Choose vector against which all angles will be measured 536 536 if len(polygon) > 1: 537 v_prev = p0  polygon[2] # Vector that leads to p0 from previous point 537 v_prev = p0  polygon[2] # Vector that leads to p0 538 # from previous point 538 539 else: 539 # FIXME (Ole): What do we do if the first point has multiple540 #candidates?540 # FIXME (Ole): What do we do if the first point has 541 # multiple candidates? 541 542 # Being the lower left corner, perhaps we can use the 542 # vector [1, 0], but I really don't know if this is completely543 #watertight.543 # vector [1, 0], but I really don't know if this is 544 # completely watertight. 544 545 v_prev = [1.0, 0.0] 545 546 … … 549 550 for pc in candidate_list: 550 551 551 552 552 vc = pcp0 # Candidate vector (from p0 to candidate pt) 553 553 … … 556 556 ac = angle(vc, v_prev) 557 557 if ac > pi: 558 # Give preference to angles on the right hand side of v_prev 559 #print 'pc = %s, changing angle from %f to %f' %(pc, ac*180/pi, (ac2*pi)*180/pi) 558 # Give preference to angles on the right hand side 559 # of v_prev 560 # print 'pc = %s, changing angle from %f to %f'\ 561 # %(pc, ac*180/pi, (ac2*pi)*180/pi) 560 562 ac = ac2*pi 561 563 562 # take the minimal angle corresponding to the rightmost vector 564 # Take the minimal angle corresponding to the 565 # rightmost vector 563 566 if ac < minimum_angle: 564 567 minimum_angle = ac … … 567 570 568 571 if verbose is True: 569 print ' Best candidate %s, angle %f' %(p1, minimum_angle*180/pi) 572 print ' Best candidate %s, angle %f'\ 573 %(p1, minimum_angle*180/pi) 570 574 571 575 else: 572 576 p1 = candidate_list[0] 573 577 578 574 579 if point_registry.has_key(tuple(p1)): 575 # We have completed the boundary polygon  yeehaa 576 break 580 # We have reached a point already visited. 581 582 if allclose(p1, polygon[0]): 583 # If it is the initial point, the polygon is complete. 584 585 if verbose is True: 586 print ' Stop criterion fulfilled at point %s' %p1 587 print polygon 588 589 # We have completed the boundary polygon  yeehaa 590 break 591 else: 592 # The point already visited is not the initial point 593 # This would be a pathological triangle, but the 594 # algorithm must be able to deal with this 595 pass 596 577 597 else: 598 # We are still finding new points on the boundary 578 599 point_registry[tuple(p1)] = len(point_registry) 579 600 
anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py
r4171 r4515 809 809 assert is_inside_polygon(p, P) 810 810 811 812 def test_boundary_polygon_IIIa(self): 813 """test_boundary_polygon_IIIa  Check pathological situation where 814 one triangle has no neighbours. This may be the case if a mesh 815 is partitioned using pymetis. 816 """ 817 818 from Numeric import zeros, Float 819 820 821 #Points 822 a = [0.0, 0.0] #0 823 b = [0.0, 0.5] #1 824 c = [0.0, 1.0] #2 825 d = [0.5, 0.0] #3 826 e = [0.5, 0.5] #4 827 f = [1.0, 0.0] #5 828 g = [1.0, 0.5] #6 829 h = [1.0, 1.0] #7 830 831 # Add pathological triangle with no neighbours to an otherwise 832 # trivial mesh 833 834 points = [a, b, c, d, e, f, g, h] 835 836 #cbe, aeb, dea, fed, ghe (pathological triangle) 837 vertices = [[2,1,4], [0,4,1], [3,4,0], [5,4,3], 838 [6,7,4]] 839 840 mesh = Mesh(points, vertices) 841 mesh.check_integrity() 842 843 P = mesh.get_boundary_polygon(verbose=False) 844 845 846 assert len(P) == 9 847 848 # Note that point e appears twice! 849 assert allclose(P, [a, d, f, e, g, h, e, c, b]) 850 851 for p in points: 852 msg = 'Point %s is not inside polygon %s'\ 853 %(p, P) 854 assert is_inside_polygon(p, P), msg 855 856 857 858 859 811 860 812 861 def test_boundary_polygon_IV(self): … … 1149 1198 # 1150 1199 if __name__ == "__main__": 1151 #suite = unittest.makeSuite(Test_Mesh,'test_ mesh_get_boundary_polygon_with_georeferencing')1200 #suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon_IIIa') 1152 1201 suite = unittest.makeSuite(Test_Mesh,'test') 1153 1202 runner = unittest.TextTestRunner()
Note: See TracChangeset
for help on using the changeset viewer.