# Changeset 5286

Ignore:
Timestamp:
May 6, 2008, 6:03:26 PM (15 years ago)
Message:

Work on ticket:175 with special consideration to mesh intersections
Also reinstated disabled Okada test which currently fails but is faster.

Location:
anuga_core/source/anuga
Files:
5 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

 r5221 # Build structure listing which trianglse belong to which nodet. # Build structure listing which trianglse belong to which node. if verbose: print 'Building inverted triangle structure' self.build_inverted_triangle_structure()
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

 r5276 from general_mesh import General_mesh from math import pi, sqrt from Numeric import array, allclose if ratio <= min_ratio: min_ratio = ratio return max_ratio,min_ratio def build_neighbour_structure(self): N = len(self) #Get x,y coordinates for all vertices for all triangles # Get x,y coordinates for all vertices for all triangles V = self.get_vertex_coordinates() #Check each triangle # Check each triangle for i in range(N): x2, y2 = V[3*i+2, :] #Check that area hasn't been compromised # Check that area hasn't been compromised area = self.areas[i] ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2 assert abs((area - ref)/area) < epsilon, msg #Check that points are arranged in counter clock-wise order # Check that points are arranged in counter clock-wise order v0 = [x1-x0, y1-y0] v1 = [x2-x1, y2-y1] assert a0 < pi and a1 < pi and a2 < pi, msg #Check that normals are orthogonal to edge vectors #Note that normal[k] lies opposite vertex k # Check that normals are orthogonal to edge vectors # Note that normal[k] lies opposite vertex k normal0 = self.normals[i, 0:2] for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]: #Normalise # Normalise l_u = sqrt(u*u + u*u) l_v = sqrt(v*v + v*v) x = (u*v + u*v)/l_u/l_v #Inner product x = (u*v + u*v)/l_u/l_v # Inner product msg = 'Normal vector (%f,%f) is not perpendicular to' %tuple(v) def _get_intersecting_segments(self, line, triangle_intersections={}): """Find edges intersected by line Input: line - list of two points forming a segmented line triangle_intersections to be updated Output: list of instances of class Triangle_intersection This method is used by the public method get_intersecting_segments(self, polyline) which also contains more documentation. """ from anuga.utilities.polygon import intersection msg = 'Line segment must contain exactly two points' assert len(line) == 2, msg # Check intersection with edge segments for all triangles # FIXME (Ole): This should be implemented in C V = self.get_vertex_coordinates() N = len(self) for i in range(N): # Get nodes and edge segments for each triangle x0, y0 = V[3*i, :] x1, y1 = V[3*i+1, :] x2, y2 = V[3*i+2, :] edge_segments = [[[x0,y0], [x1, y1]], [[x1,y1], [x2, y2]], [[x2,y2], [x0, y0]]] # Find segments that are intersected by line intersections = [] for edge in edge_segments: status, value = intersection(line, edge) if status == 1: # Normal intersection of one edge # Exclude singular intersections with vertices if not(allclose(value, edge) or\ allclose(value, edge)): intersections.append(value) if status == 2: # Edge is sharing a segment with line # Add identified line segment for j in range(value.shape): intersections.append(value[j,:]) msg = 'There can be only two or no intersections' assert len(intersections) in [0,2], msg if len(intersections) == 2: # Calculate attributes for this segment # Origin of intersecting line to be used for direction xi0 = line eta0 = line # End points of intersecting segment x0, y0 = intersections x1, y1 = intersections # Determine which end point is closer to the origin of the line # This is necessary for determining the direction of # the line and the normals # Distances from line origin to the two intersections z0 = array([x0 - xi0, y0 - eta0]) z1 = array([x1 - xi0, y1 - eta0]) d0 = sqrt(sum(z0**2)) d1 = sqrt(sum(z1**2)) if d1 < d0: # Swap xi, eta = x0, y0 x0, y0 = x1, y1 x1, y1 = xi, eta # (x0,y0) is now the origin of the intersecting segment # Normal direction (right hand side relative to direction of line) vector = array([x1 - x0, y1 - y0]) # Segment vector length = sqrt(sum(vector**2))      # Segment length normal = array([vector, -vector])/length segment = ((x0,y0), (x1, y1)) T = Triangle_intersection(segment=segment, normal=normal, length=length, triangle_id=i) # Add segment unless it was done earlier if not triangle_intersections.has_key(segment): triangle_intersections[segment] = T #return triangle_intersections def get_intersecting_segments(self, polyline): """Find edges intersected by polyline Output: dictionary where keys are triangle ids for those elements that are intersected values are a list of segments each represented as a triplet: length of the intersecting segment right hand normal midpoint of intersecting segment The polyline may break inside any triangle causing multiple segments per triabgle. list of instances of class Triangle_intersection The polyline may break inside any triangle causing multiple segments per triangle - consequently the same triangle may appear in several entries. If a polyline segment coincides with a triangle edge, the the entire shared segment will be used. Onle one of the triangles thus intersected will be used and that is the first one encoutered. intersections with single vertices are ignored. Resulting segments are unsorted """ pass msg = 'Polyline must contain at least two points' assert len(polyline) >= 2, msg # For all segments in polyline triangle_intersections = {} for i, point0 in enumerate(polyline[:-1]): point1 = polyline[i+1] line = [point0, point1] self._get_intersecting_segments(line, triangle_intersections) return triangle_intersections.values() return [] class Triangle_intersection: """Store information about line segments intersecting a triangle Attributes are segment: Line segment intersecting triangle [[x0,y0], [x1, y1]] normal: [a,b] right hand normal to segment length: Length of intersecting segment triangle_id: id (in mesh) of triangle being intersected """ def __init__(self, segment=None, normal=None, length=None, triangle_id=None): self.segment = segment self.normal = normal self.length = length self.triangle_id = triangle_id def __repr__(self): s = 'Triangle_intersection(' s += 'segment=%s, normal=%s, length=%s, triangle_id=%s)'\ %(self.segment, self.normal, self.length, self.triangle_id) return s
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

 r5276 def NOtest_get_intersecting_segments(self): def test_get_intersecting_segments1(self): """test_get_intersecting_segments(self): Very simple test """ # Build test mesh # Create basic mesh # 9 points at (0,0), (0, 1), ..., (2,2) # 8 triangles enumerated from left bottom to right top. points, vertices, boundary = rectangular(2, 2, 2, 2) mesh = Mesh(points, vertices, boundary) # Very simple horizontal line intersecting # for y_line in [0.1, 0.2, 0.314159, 0.41, 0.6, 0.99, 1.01, 1.5, 1.77, 1.9]: if y_line < 1: ceiling = 1 floor = 0 intersected_triangles = [0,1,4,5] elif y_line > 1: ceiling = 2 floor = 1 intersected_triangles = [2,3,6,7] else: raise Exception, 'this test is not for parallel lines' line = [[-1,y_line], [3,y_line]] L = mesh.get_intersecting_segments(line) assert len(L) == 4 # Check all normals point straight down etc for x in L: if x.triangle_id % 2 == 0: assert allclose(x.length, ceiling-y_line) else: assert allclose(x.length, y_line-floor) assert allclose(x.normal, [0,-1]) assert allclose(x.segment, x.segment + x.length) assert allclose(x.segment, y_line) assert allclose(x.segment, y_line) assert x.triangle_id in intersected_triangles def test_get_intersecting_segments_coinciding(self): """test_get_intersecting_segments_coinciding(self): Test that lines coinciding with triangle edges work. """ # Build test mesh # Create basic mesh # 9 points at (0,0), (0, 1), ..., (2,2) # 8 triangles enumerated from left bottom to right top. points, vertices, boundary = rectangular(2, 2, 2, 2) mesh = Mesh(points, vertices, boundary) intersected_triangles = [1,5] # Very simple horizontal line intersecting # y_line = 1.0 line = [[-1,y_line], [3,y_line]] L = mesh.get_intersecting_segments(line) msg = 'Only two triangles should be returned' assert len(L) == 2, msg # Check all normals point straight down for i, x in enumerate(L): assert allclose(x.length, 1.0) assert allclose(x.normal, [0,-1]) assert allclose(x.segment, x.segment + x.length) assert allclose(x.segment, y_line) assert allclose(x.segment, y_line) assert x.triangle_id in intersected_triangles def test_get_intersecting_segments2(self): """test_get_intersecting_segments(self): More complex mesh """ # Build test mesh (from bounding polygon tests pass if __name__ == "__main__": #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles') #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_coinciding') suite = unittest.makeSuite(Test_Mesh,'test') runner = unittest.TextTestRunner()