- Timestamp:
- Feb 27, 2009, 11:54:09 AM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r6304 r6428 13 13 14 14 import numpy as num 15 15 16 16 17 17 class Mesh(General_mesh): … … 91 91 verbose=verbose) 92 92 93 if verbose: print 'Initialising mesh' 93 if verbose: print 'Initialising mesh' 94 94 95 95 N = len(self) #Number_of_triangles … … 111 111 112 112 #Initialise each triangle 113 if verbose: print 'Mesh: Computing centroids and radii' 113 if verbose: print 'Mesh: Computing centroids and radii' 114 114 for i in range(N): 115 115 if verbose and i % ((N+10)/10) == 0: print '(%d/%d)' %(i, N) … … 117 117 x0, y0 = V[3*i, :] 118 118 x1, y1 = V[3*i+1, :] 119 x2, y2 = V[3*i+2, :] 119 x2, y2 = V[3*i+2, :] 120 120 121 121 #x0 = V[i, 0]; y0 = V[i, 1] … … 166 166 167 167 #Build neighbour structure 168 if verbose: print 'Mesh: Building neigbour structure' 168 if verbose: print 'Mesh: Building neigbour structure' 169 169 self.build_neighbour_structure() 170 170 … … 178 178 179 179 #Build tagged element dictionary mapping (tag) to array of elements 180 if verbose: print 'Mesh: Building tagged elements dictionary' 180 if verbose: print 'Mesh: Building tagged elements dictionary' 181 181 self.build_tagged_elements_dictionary(tagged_elements) 182 182 … … 184 184 self.lone_vertices = [] 185 185 #Check that all vertices have been registered 186 for node, count in enumerate(self.number_of_triangles_per_node): 186 for node, count in enumerate(self.number_of_triangles_per_node): 187 187 #msg = 'Node %d does not belong to an element.' %node 188 188 #assert count > 0, msg 189 189 if count == 0: 190 190 self.lone_vertices.append(node) 191 191 192 192 #Update boundary indices FIXME: OBSOLETE 193 193 #self.build_boundary_structure() 194 194 195 195 #FIXME check integrity? 196 if verbose: print 'Mesh: Done' 196 if verbose: print 'Mesh: Done' 197 197 198 198 def __repr__(self): … … 400 400 #print "tagged_elements", tagged_elements 401 401 self.tagged_elements = tagged_elements 402 402 403 403 def get_tagged_elements(self): 404 404 return self.tagged_elements … … 459 459 460 460 Using the mesh boundary, derive a bounding polygon for this mesh. 461 If multiple vertex values are present (vertices stored uniquely), 462 461 If multiple vertex values are present (vertices stored uniquely), 462 the algorithm will select the path that contains the entire mesh. 463 463 464 464 All points are in absolute UTM coordinates 465 465 """ 466 467 from anuga.utilities.numerical_tools import angle, ensure_numeric 468 466 467 from anuga.utilities.numerical_tools import angle, ensure_numeric 469 468 470 469 # Get mesh extent 471 470 xmin, xmax, ymin, ymax = self.get_extent(absolute=True) 472 471 pmin = ensure_numeric([xmin, ymin]) 473 pmax = ensure_numeric([xmax, ymax]) 474 472 pmax = ensure_numeric([xmax, ymax]) 475 473 476 474 # Assemble dictionary of boundary segments and choose starting point … … 478 476 inverse_segments = {} 479 477 p0 = None 480 mindist = num.sqrt(num.sum((pmax-pmin)**2)) # Start value across entire mesh 478 479 # Start value across entire mesh 480 mindist = num.sqrt(num.sum((pmax-pmin)**2)) 481 481 for i, edge_id in self.boundary.keys(): 482 482 # Find vertex ids for boundary segment … … 485 485 if edge_id == 2: a = 0; b = 1 486 486 487 A = self.get_vertex_coordinate(i, a, absolute=True) # Start 488 B = self.get_vertex_coordinate(i, b, absolute=True) # End 489 487 A = self.get_vertex_coordinate(i, a, absolute=True) # Start 488 B = self.get_vertex_coordinate(i, b, absolute=True) # End 490 489 491 490 # Take the point closest to pmin as starting point … … 503 502 p0 = B 504 503 505 506 504 # Sanity check 507 505 if p0 is None: … … 511 509 # Register potential paths from A to B 512 510 if not segments.has_key(tuple(A)): 513 segments[tuple(A)] = [] # Empty list for candidate points 514 515 segments[tuple(A)].append(B) 516 511 segments[tuple(A)] = [] # Empty list for candidate points 512 513 segments[tuple(A)].append(B) 517 514 518 515 # Start with smallest point and follow boundary (counter clock wise) 519 516 polygon = [list(p0)]# Storage for final boundary polygon 520 point_registry = {} # Keep track of storage to avoid multiple runs 521 # around boundary. This will only be the case if 522 517 point_registry = {} # Keep track of storage to avoid multiple runs 518 # around boundary. This will only be the case if 519 # there are more than one candidate. 523 520 # FIXME (Ole): Perhaps we can do away with polygon 524 521 # and use only point_registry to save space. 525 522 526 point_registry[tuple(p0)] = 0 527 523 point_registry[tuple(p0)] = 0 524 528 525 while len(point_registry) < len(self.boundary): 529 530 526 candidate_list = segments[tuple(p0)] 531 527 if len(candidate_list) > 1: 532 # Multiple points detected (this will be the case for meshes 533 # with duplicate points as those used for discontinuous 534 535 # Take the candidate that is furthest to the clockwise 536 537 538 539 528 # Multiple points detected (this will be the case for meshes 529 # with duplicate points as those used for discontinuous 530 # triangles with vertices stored uniquely). 531 # Take the candidate that is furthest to the clockwise 532 # direction, as that will follow the boundary. 533 # 534 # This will also be the case for pathological triangles 535 # that have no neighbours. 540 536 541 537 if verbose: 542 print 'Point %s has multiple candidates: %s'\543 %(str(p0), candidate_list)538 print ('Point %s has multiple candidates: %s' 539 % (str(p0), candidate_list)) 544 540 545 541 # Check that previous are not in candidate list … … 548 544 549 545 # Choose vector against which all angles will be measured 550 if len(polygon) > 1: 551 v_prev = p0 - polygon[-2] # Vector that leads to p0552 546 if len(polygon) > 1: 547 v_prev = p0 - polygon[-2] # Vector that leads to p0 548 # from previous point 553 549 else: 554 # FIXME (Ole): What do we do if the first point has 555 550 # FIXME (Ole): What do we do if the first point has 551 # multiple candidates? 556 552 # Being the lower left corner, perhaps we can use the 557 # vector [1, 0], but I really don't know if this is 558 553 # vector [1, 0], but I really don't know if this is 554 # completely watertight. 559 555 v_prev = [1.0, 0.0] 560 561 562 # Choose candidate with minimum angle 556 557 # Choose candidate with minimum angle 563 558 minimum_angle = 2*pi 564 559 for pc in candidate_list: 565 566 vc = pc-p0 # Candidate vector (from p0 to candidate pt) 567 560 vc = pc-p0 # Candidate vector (from p0 to candidate pt) 561 568 562 # Angle between each candidate and the previous vector 569 563 # in [-pi, pi] 570 564 ac = angle(vc, v_prev) 571 565 if ac > pi: 572 # Give preference to angles on the right hand side 573 # of v_prev 566 # Give preference to angles on the right hand side 567 # of v_prev 574 568 # print 'pc = %s, changing angle from %f to %f'\ 575 569 # %(pc, ac*180/pi, (ac-2*pi)*180/pi) 576 570 ac = ac-2*pi 577 571 578 # Take the minimal angle corresponding to the 579 572 # Take the minimal angle corresponding to the 573 # rightmost vector 580 574 if ac < minimum_angle: 581 575 minimum_angle = ac 582 p1 = pc # Best candidate 583 576 p1 = pc # Best candidate 584 577 585 578 if verbose is True: 586 579 print ' Best candidate %s, angle %f'\ 587 580 %(p1, minimum_angle*180/pi) 588 589 581 else: 590 582 p1 = candidate_list[0] 591 583 592 593 584 if point_registry.has_key(tuple(p1)): 594 # We have reached a point already visited. 595 596 if num.allclose(p1, polygon[0]): 597 # If it is the initial point, the polygon is complete. 598 585 # We have reached a point already visited. 586 if num.allclose(p1, polygon[0]): 587 # If it is the initial point, the polygon is complete. 599 588 if verbose is True: 600 601 print polygon 602 589 print ' Stop criterion fulfilled at point %s' %p1 590 print polygon 591 603 592 # We have completed the boundary polygon - yeehaa 604 605 else: 606 607 # This would be a pathological triangle, but the 608 609 610 593 break 594 else: 595 # The point already visited is not the initial point 596 # This would be a pathological triangle, but the 597 # algorithm must be able to deal with this 598 pass 599 611 600 else: 612 601 # We are still finding new points on the boundary 613 602 point_registry[tuple(p1)] = len(point_registry) 614 615 polygon.append(list(p1)) # De-numeric each point :-)603 604 polygon.append(list(p1)) # De-numeric each point :-) 616 605 p0 = p1 617 606 618 619 607 return polygon 620 621 608 622 609 def check_integrity(self): … … 641 628 x1, y1 = V[3*i+1, :] 642 629 x2, y2 = V[3*i+2, :] 643 630 644 631 # Check that area hasn't been compromised 645 632 area = self.areas[i] … … 678 665 679 666 x = (u[0]*v[0] + u[1]*v[1])/l_u # Inner product 680 667 681 668 msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v) 682 669 msg += ' edge (%f,%f) in triangle %d.' %(tuple(u) + (i,)) 683 msg += ' Inner product is %e.' %x 670 msg += ' Inner product is %e.' %x 684 671 assert x < epsilon, msg 685 672 … … 690 677 for i in range(N): 691 678 # For each triangle 692 679 693 680 for k, neighbour_id in enumerate(self.neighbours[i,:]): 694 681 … … 752 739 # Node is lone - i.e. not part of the mesh 753 740 continue 754 741 755 742 k += 1 756 743 757 744 volume_id = index / 3 758 745 vertex_id = index % 3 759 746 760 747 msg = 'Triangle %d, vertex %d points to %d. Should have been %d'\ 761 748 %(volume_id, vertex_id, self.triangles[volume_id, vertex_id], current_node) 762 749 assert self.triangles[volume_id, vertex_id] == current_node, msg 763 750 764 751 if self.number_of_triangles_per_node[current_node] == k: 765 752 # Move on to next node … … 788 775 if not self.geo_reference.is_absolute(): 789 776 V = self.geo_reference.get_absolute(V) 790 777 791 778 return V 792 779 793 780 794 781 def get_radii(self): 795 782 """Return all radii. 796 783 Return radius of inscribed cirle for all triangles 797 784 """ 798 return self.radii 785 return self.radii 799 786 800 787 … … 826 813 str += ' Areas [m^2]:\n' 827 814 str += ' A in [%f, %f]\n' %(min(areas), max(areas)) 828 str += ' number of distinct areas: %d\n' %(len(areas)) 815 str += ' number of distinct areas: %d\n' %(len(areas)) 829 816 str += ' Histogram:\n' 830 817 … … 833 820 lo = hi 834 821 if i+1 < len(bins): 835 #Open upper interval 822 #Open upper interval 836 823 hi = bins[i+1] 837 str += ' [%f, %f[: %d\n' %(lo, hi, count) 824 str += ' [%f, %f[: %d\n' %(lo, hi, count) 838 825 else: 839 826 #Closed upper interval … … 849 836 k = 0 850 837 lower = min(areas) 851 for i, a in enumerate(areas): 852 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas 838 for i, a in enumerate(areas): 839 if i % (N/10) == 0 and i != 0: #For every 10% of the sorted areas 853 840 str += ' %d triangles in [%f, %f]\n' %(i-k, lower, a) 854 841 lower = a 855 842 k = i 856 843 857 844 str += ' %d triangles in [%f, %f]\n'\ 858 %(N-k, lower, max(areas)) 859 860 845 %(N-k, lower, max(areas)) 846 847 861 848 str += ' Boundary:\n' 862 849 str += ' Number of boundary segments == %d\n' %(len(self.boundary)) 863 str += ' Boundary tags == %s\n' %self.get_boundary_tags() 850 str += ' Boundary tags == %s\n' %self.get_boundary_tags() 864 851 str += '------------------------------------------------\n' 865 852 866 853 867 854 return str 868 855 869 856 870 857 def get_triangle_containing_point(self, point): … … 874 861 875 862 """ 876 863 877 864 # FIXME(Ole): This function is currently brute force 878 865 # because I needed it for diagnostics. … … 884 871 polygon = self.get_boundary_polygon() 885 872 #print polygon, point 886 873 887 874 if is_outside_polygon(point, polygon): 888 875 msg = 'Point %s is outside mesh' %str(point) … … 899 886 if is_inside_polygon(point, poly, closed=True): 900 887 return i 901 888 902 889 return 903 890 … … 930 917 Resulting segments are unsorted 931 918 """ 932 919 933 920 V = self.get_vertex_coordinates() 934 921 N = len(self) 935 922 936 923 # Adjust polyline to mesh spatial origin 937 924 polyline = self.geo_reference.get_relative(polyline) … … 939 926 if use_cache is True: 940 927 segments = cache(get_intersecting_segments, 941 (V, N, polyline), 928 (V, N, polyline), 942 929 {'verbose': verbose}, 943 930 verbose=verbose) 944 else: 931 else: 945 932 segments = get_intersecting_segments(V, N, polyline, 946 933 verbose=verbose) 947 934 948 935 949 936 return segments 950 951 937 938 952 939 953 940 def get_triangle_neighbours(self, tri_id): … … 956 943 957 944 Negative returned triangle id's represent a boundary as a neighbour. 958 945 959 946 If the given triangle id is bad, return an empty list. 960 947 """ … … 964 951 except IndexError: 965 952 return [] 966 953 967 954 968 955 def get_interpolation_object(self): 969 956 """Get object I that will allow linear interpolation using this mesh 970 971 This is a time consuming process but it needs only to be 957 958 This is a time consuming process but it needs only to be 972 959 once for the mesh. 973 974 Interpolation can then be done using 975 976 result = I.interpolate_block(vertex_values, interpolation_points) 977 960 961 Interpolation can then be done using 962 963 result = I.interpolate_block(vertex_values, interpolation_points) 964 978 965 where vertex values have been obtained from a quantity using 979 966 vertex_values, triangles = self.get_vertex_values() … … 984 971 else: 985 972 from anuga.fit_interpolate.interpolate import Interpolate 986 987 # Get discontinuous mesh - this will match internal 973 974 # Get discontinuous mesh - this will match internal 988 975 # representation of vertex values 989 976 triangles = self.get_disconnected_triangles() … … 992 979 I = Interpolate(vertex_coordinates, triangles) 993 980 self.interpolation_object = I 994 995 return I 996 981 982 return I 983 997 984 998 985 class Triangle_intersection: 999 986 """Store information about line segments intersecting a triangle 1000 987 1001 988 Attributes are 1002 989 … … 1014 1001 length=None, 1015 1002 triangle_id=None): 1016 self.segment = segment 1003 self.segment = segment 1017 1004 self.normal = normal 1018 1005 self.length = length 1019 1006 self.triangle_id = triangle_id 1020 1007 1021 1008 1022 1009 def __repr__(self): … … 1027 1014 self.length, 1028 1015 self.triangle_id) 1029 1016 1030 1017 return s 1031 1018 1032 1019 1033 1020 1034 1021 def _get_intersecting_segments(V, N, line, … … 1051 1038 from anuga.utilities.polygon import intersection 1052 1039 from anuga.utilities.polygon import is_inside_polygon 1053 1040 1054 1041 msg = 'Line segment must contain exactly two points' 1055 1042 assert len(line) == 2, msg … … 1060 1047 eta0 = line[0][1] 1061 1048 1062 1049 1063 1050 # Check intersection with edge segments for all triangles 1064 1051 # FIXME (Ole): This should be implemented in C … … 1069 1056 x1, y1 = V[3*i+1, :] 1070 1057 x2, y2 = V[3*i+2, :] 1071 1058 1072 1059 1073 1060 edge_segments = [[[x0,y0], [x1, y1]], … … 1076 1063 1077 1064 # Find segments that are intersected by line 1078 1065 1079 1066 intersections = {} # Use dictionary to record points only once 1080 1067 for edge in edge_segments: … … 1082 1069 status, value = intersection(line, edge) 1083 1070 #if value is not None: print 'Triangle %d, Got' %i, status, value 1084 1071 1085 1072 if status == 1: 1086 1073 # Normal intersection of one edge or vertex 1087 intersections[tuple(value)] = i 1074 intersections[tuple(value)] = i 1088 1075 1089 1076 # Exclude singular intersections with vertices … … 1101 1088 # along this edge, it will be recorded here. 1102 1089 intersections[tuple(value[0,:])] = i 1103 intersections[tuple(value[1,:])] = i 1104 1105 1090 intersections[tuple(value[1,:])] = i 1091 1092 1106 1093 if len(intersections) == 1: 1107 1094 # Check if either line end point lies fully within this triangle … … 1113 1100 intersections[tuple(line[1])] = i 1114 1101 elif is_inside_polygon(line[0], poly, closed=False): 1115 intersections[tuple(line[0])] = i 1102 intersections[tuple(line[0])] = i 1116 1103 else: 1117 # Ignore situations where one vertex is touch, for instance 1104 # Ignore situations where one vertex is touch, for instance 1118 1105 continue 1119 1106 … … 1141 1128 z0 = num.array([x0 - xi0, y0 - eta0], num.float) 1142 1129 z1 = num.array([x1 - xi0, y1 - eta0], num.float) 1143 d0 = num.sqrt(num.sum(z0**2)) 1130 d0 = num.sqrt(num.sum(z0**2)) 1144 1131 d1 = num.sqrt(num.sum(z1**2)) 1145 1132 1146 1133 if d1 < d0: 1147 1134 # Swap … … 1151 1138 1152 1139 # (x0,y0) is now the origin of the intersecting segment 1153 1140 1154 1141 1155 1142 # Normal direction: … … 1160 1147 1161 1148 1162 segment = ((x0,y0), (x1, y1)) 1149 segment = ((x0,y0), (x1, y1)) 1163 1150 T = Triangle_intersection(segment=segment, 1164 1151 normal=normal, … … 1168 1155 1169 1156 # Add segment unless it was done earlier 1170 if not triangle_intersections.has_key(segment): 1157 if not triangle_intersections.has_key(segment): 1171 1158 triangle_intersections[segment] = T 1172 1159 1173 1160 1174 # Return segments as a list 1161 # Return segments as a list 1175 1162 return triangle_intersections.values() 1176 1163 1177 1164 1178 1165 def get_intersecting_segments(V, N, polyline, 1179 verbose=False): 1166 verbose=False): 1180 1167 """Internal function to find edges intersected by Polyline 1181 1168 1182 1169 Input: 1183 1170 V: Vertex coordinates as obtained by mesh.get_vertex_coordinates() 1184 1171 N: Number of triangles in mesh 1185 polyline - list of points forming a segmented line 1172 polyline - list of points forming a segmented line 1186 1173 verbose 1187 1174 Output: … … 1190 1177 This method is used by the public method 1191 1178 get_intersecting_segments(self, polyline) which also contains 1192 more documentation. 1179 more documentation. 1193 1180 """ 1194 1181 1195 1182 msg = 'Polyline must contain at least two points' 1196 1183 assert len(polyline) >= 2, msg 1197 1198 1184 1185 1199 1186 # For all segments in polyline 1200 1187 triangle_intersections = [] … … 1206 1193 print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0[0], point0[1], 1207 1194 point1[0], point1[1]) 1208 1195 1209 1196 line = [point0, point1] 1210 1197 triangle_intersections += _get_intersecting_segments(V, N, line, … … 1214 1201 msg = 'No segments found' 1215 1202 assert len(triangle_intersections) > 0, msg 1216 1217 1203 1204 1218 1205 return triangle_intersections 1219 1220 1221 1222 1223 1206 1207 1208 1209 1210 1224 1211 def segment_midpoints(segments): 1225 1212 """Calculate midpoints of all segments 1226 1213 1227 1214 Inputs: 1228 1215 segments: List of instances of class Segment 1229 1216 1230 1217 Ouputs: 1231 1218 midpoints: List of points 1232 1219 """ 1233 1220 1234 1221 midpoints = [] 1235 1222 msg = 'Elements of input list to segment_midpoints must be of class Triangle_intersection' 1236 1223 for segment in segments: 1237 1224 assert isinstance(segment, Triangle_intersection), msg 1238 1239 midpoint = num.sum(num.array(segment.segment, num.float) )/21225 1226 midpoint = num.sum(num.array(segment.segment, num.float), axis=0)/2 1240 1227 midpoints.append(midpoint) 1241 1228 1242 1229 return midpoints 1243 1244 1245 1230 1231 1232
Note: See TracChangeset
for help on using the changeset viewer.