 Timestamp:
 May 6, 2008, 6:03:26 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py
r5276 r5286 10 10 from general_mesh import General_mesh 11 11 from math import pi, sqrt 12 from Numeric import array, allclose 12 13 13 14 … … 230 231 if ratio <= min_ratio: min_ratio = ratio 231 232 return max_ratio,min_ratio 233 234 232 235 233 236 def build_neighbour_structure(self): … … 628 631 629 632 N = len(self) 630 #Get x,y coordinates for all vertices for all triangles 633 634 # Get x,y coordinates for all vertices for all triangles 631 635 V = self.get_vertex_coordinates() 632 #Check each triangle 636 637 # Check each triangle 633 638 for i in range(N): 634 639 … … 637 642 x2, y2 = V[3*i+2, :] 638 643 639 # Check that area hasn't been compromised644 # Check that area hasn't been compromised 640 645 area = self.areas[i] 641 646 ref = abs((x1*y0x0*y1)+(x2*y1x1*y2)+(x0*y2x2*y0))/2 … … 644 649 assert abs((area  ref)/area) < epsilon, msg 645 650 646 # Check that points are arranged in counter clockwise order651 # Check that points are arranged in counter clockwise order 647 652 v0 = [x1x0, y1y0] 648 653 v1 = [x2x1, y2y1] … … 656 661 assert a0 < pi and a1 < pi and a2 < pi, msg 657 662 658 # Check that normals are orthogonal to edge vectors659 # Note that normal[k] lies opposite vertex k663 # Check that normals are orthogonal to edge vectors 664 # Note that normal[k] lies opposite vertex k 660 665 661 666 normal0 = self.normals[i, 0:2] … … 665 670 for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]: 666 671 667 # Normalise672 # Normalise 668 673 l_u = sqrt(u[0]*u[0] + u[1]*u[1]) 669 674 l_v = sqrt(v[0]*v[0] + v[1]*v[1]) 670 675 671 x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v # Inner product676 x = (u[0]*v[0] + u[1]*v[1])/l_u/l_v # Inner product 672 677 673 678 msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v) … … 896 901 897 902 903 def _get_intersecting_segments(self, line, triangle_intersections={}): 904 """Find edges intersected by line 905 906 Input: 907 line  list of two points forming a segmented line 908 triangle_intersections to be updated 909 Output: 910 list of instances of class Triangle_intersection 911 912 This method is used by the public method 913 get_intersecting_segments(self, polyline) which also contains 914 more documentation. 915 """ 916 917 from anuga.utilities.polygon import intersection 918 919 msg = 'Line segment must contain exactly two points' 920 assert len(line) == 2, msg 921 922 923 # Check intersection with edge segments for all triangles 924 # FIXME (Ole): This should be implemented in C 925 V = self.get_vertex_coordinates() 926 N = len(self) 927 928 for i in range(N): 929 # Get nodes and edge segments for each triangle 930 x0, y0 = V[3*i, :] 931 x1, y1 = V[3*i+1, :] 932 x2, y2 = V[3*i+2, :] 933 934 935 edge_segments = [[[x0,y0], [x1, y1]], 936 [[x1,y1], [x2, y2]], 937 [[x2,y2], [x0, y0]]] 938 939 # Find segments that are intersected by line 940 intersections = [] 941 for edge in edge_segments: 942 943 status, value = intersection(line, edge) 944 if status == 1: 945 # Normal intersection of one edge 946 947 # Exclude singular intersections with vertices 948 if not(allclose(value, edge[0]) or\ 949 allclose(value, edge[1])): 950 intersections.append(value) 951 952 if status == 2: 953 # Edge is sharing a segment with line 954 955 # Add identified line segment 956 for j in range(value.shape[0]): 957 intersections.append(value[j,:]) 958 959 960 msg = 'There can be only two or no intersections' 961 assert len(intersections) in [0,2], msg 962 963 964 if len(intersections) == 2: 965 966 # Calculate attributes for this segment 967 968 # Origin of intersecting line to be used for direction 969 xi0 = line[0][0] 970 eta0 = line[0][1] 971 972 # End points of intersecting segment 973 x0, y0 = intersections[0] 974 x1, y1 = intersections[1] 975 976 977 # Determine which end point is closer to the origin of the line 978 # This is necessary for determining the direction of 979 # the line and the normals 980 981 # Distances from line origin to the two intersections 982 z0 = array([x0  xi0, y0  eta0]) 983 z1 = array([x1  xi0, y1  eta0]) 984 d0 = sqrt(sum(z0**2)) 985 d1 = sqrt(sum(z1**2)) 986 987 if d1 < d0: 988 # Swap 989 xi, eta = x0, y0 990 x0, y0 = x1, y1 991 x1, y1 = xi, eta 992 993 # (x0,y0) is now the origin of the intersecting segment 994 995 996 # Normal direction (right hand side relative to direction of line) 997 vector = array([x1  x0, y1  y0]) # Segment vector 998 length = sqrt(sum(vector**2)) # Segment length 999 normal = array([vector[1], vector[0]])/length 1000 1001 1002 segment = ((x0,y0), (x1, y1)) 1003 T = Triangle_intersection(segment=segment, 1004 normal=normal, 1005 length=length, 1006 triangle_id=i) 1007 1008 1009 # Add segment unless it was done earlier 1010 if not triangle_intersections.has_key(segment): 1011 triangle_intersections[segment] = T 1012 1013 1014 #return triangle_intersections 1015 1016 1017 898 1018 def get_intersecting_segments(self, polyline): 899 1019 """Find edges intersected by polyline … … 903 1023 904 1024 Output: 905 dictionary where 906 keys are triangle ids for those elements that are intersected 907 values are a list of segments each represented as a triplet: 908 length of the intersecting segment 909 right hand normal 910 midpoint of intersecting segment 911 912 The polyline may break inside any triangle causing multiple segments per triabgle. 1025 list of instances of class Triangle_intersection 1026 1027 The polyline may break inside any triangle causing multiple 1028 segments per triangle  consequently the same triangle may 1029 appear in several entries. 1030 1031 If a polyline segment coincides with a triangle edge, 1032 the the entire shared segment will be used. 1033 1034 Onle one of the triangles thus intersected will be used and that 1035 is the first one encoutered. 1036 1037 intersections with single vertices are ignored. 1038 1039 Resulting segments are unsorted 913 1040 """ 914 1041 915 pass 1042 msg = 'Polyline must contain at least two points' 1043 assert len(polyline) >= 2, msg 1044 1045 # For all segments in polyline 1046 triangle_intersections = {} 1047 for i, point0 in enumerate(polyline[:1]): 1048 point1 = polyline[i+1] 1049 1050 line = [point0, point1] 1051 1052 self._get_intersecting_segments(line, triangle_intersections) 1053 1054 1055 return triangle_intersections.values() 1056 916 1057 917 1058 … … 931 1072 return [] 932 1073 1074 1075 1076 class Triangle_intersection: 1077 """Store information about line segments intersecting a triangle 1078 1079 Attributes are 1080 1081 segment: Line segment intersecting triangle [[x0,y0], [x1, y1]] 1082 normal: [a,b] right hand normal to segment 1083 length: Length of intersecting segment 1084 triangle_id: id (in mesh) of triangle being intersected 1085 1086 """ 1087 1088 1089 def __init__(self, 1090 segment=None, 1091 normal=None, 1092 length=None, 1093 triangle_id=None): 1094 self.segment = segment 1095 self.normal = normal 1096 self.length = length 1097 self.triangle_id = triangle_id 1098 1099 1100 def __repr__(self): 1101 s = 'Triangle_intersection(' 1102 s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\ 1103 %(self.segment, 1104 self.normal, 1105 self.length, 1106 self.triangle_id) 1107 1108 return s
Note: See TracChangeset
for help on using the changeset viewer.