Changeset 5287 for anuga_core/source/anuga

Timestamp:

May 7, 2008, 12:30:35 PM (17 years ago)

Author:

ole

Message:

More work on mesh intersections and much more testing - this includesd
pathological cases where lines intersect vertices, edges and where a
polyline breaks inside a triangle.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• neighbour_mesh.py (modified) (10 diffs)
• test_neighbour_mesh.py (modified) (6 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -901 +901 @@
 
 
-    def _get_intersecting_segments(self, line, triangle_intersections={}):
+    def _get_intersecting_segments(self, line):
       """Find edges intersected by line
 
       Input:
-          line - list of two points forming a segmented line
-          triangle_intersections to be updated
+          line - list of two points forming a segmented line  
       Output:
           list of instances of class Triangle_intersection
@@ -916 +915 @@
 
       from anuga.utilities.polygon import intersection
+      from anuga.utilities.polygon import is_inside_polygon
       
       msg = 'Line segment must contain exactly two points'
       assert len(line) == 2, msg
-      
+
+      # Origin of intersecting line to be used for
+      # establishing direction
+      xi0 = line[0][0]
+      eta0 = line[0][1]
+
       
       # Check intersection with edge segments for all triangles
@@ -925 +930 @@
       V = self.get_vertex_coordinates()
       N = len(self)
-
+      triangle_intersections={} # Keep track of segments already done
       for i in range(N):
           # Get nodes and edge segments for each triangle
@@ -938 +943 @@
 
           # Find segments that are intersected by line
-          intersections = []
+          
+          intersections = {} # Use dictionary to record points only once
           for edge in edge_segments:
 
               status, value = intersection(line, edge)
               if status == 1:
-                  # Normal intersection of one edge
+                  # Normal intersection of one edge or vertex
+                  intersections[tuple(value)] = i                  
 
                   # Exclude singular intersections with vertices
-                  if not(allclose(value, edge[0]) or\
-                         allclose(value, edge[1])):
-                      intersections.append(value)
+                  #if not(allclose(value, edge[0]) or\
+                  #       allclose(value, edge[1])):
+                  #    intersections.append(value)
 
               if status == 2:
                   # Edge is sharing a segment with line
 
-                  # Add identified line segment
-                  for j in range(value.shape[0]):
-                      intersections.append(value[j,:])
+                  # This is currently covered by the two
+                  # vertices that would have been picked up
+                  # under status == 1
+                  pass
+                  
+
+
+          if len(intersections) == 1:
+              # Check if either line end point lies fully within this triangle
+              # If this is the case accept that as one end of the intersecting
+              # segment
+
+              poly = V[3*i:3*i+3]
+              if is_inside_polygon(line[1], poly, closed=False):
+                  intersections[tuple(line[1])] = i
+              elif is_inside_polygon(line[0], poly, closed=False):
+                  intersections[tuple(line[0])] = i         
+              else:
+                  # Ignore situations where one vertex is touch, for instance                   
+                  continue
 
 
@@ -966 +990 @@
               # Calculate attributes for this segment
 
-              # Origin of intersecting line to be used for direction
-              xi0 = line[0][0]
-              eta0 = line[0][1]
 
               # End points of intersecting segment
-              x0, y0 = intersections[0]
-              x1, y1 = intersections[1]
+              points = intersections.keys()
+              x0, y0 = points[0]
+              x1, y1 = points[1]
 
 
@@ -994 +1016 @@
                   
 
-              # Normal direction (right hand side relative to direction of line)
+              # Normal direction:
+              # Right hand side relative to line direction
               vector = array([x1 - x0, y1 - y0]) # Segment vector
               length = sqrt(sum(vector**2))      # Segment length
@@ -1012 +1035 @@
 
 
-      #return triangle_intersections
+      # Return segments as a list            
+      return triangle_intersections.values()
 
 
@@ -1031 +1055 @@
       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.
+      Intersections with single vertices are ignored.
 
       Resulting segments are unsorted
@@ -1044 +1067 @@
 
       # For all segments in polyline
-      triangle_intersections = {}
+      triangle_intersections = []
       for i, point0 in enumerate(polyline[:-1]):
           point1 = polyline[i+1]
@@ -1050 +1073 @@
           line = [point0, point1]
 
-          self._get_intersecting_segments(line, triangle_intersections)
-
-
-      return triangle_intersections.values()
+          triangle_intersections += self._get_intersecting_segments(line)
+
+
+      return triangle_intersections
   
 

anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

@@ -1251 +1251 @@
         """test_get_intersecting_segments(self):
 
-        Very simple test
+        Very simple test (horizontal lines)
         """
 
@@ -1287 +1287 @@
 
             # Check all normals point straight down etc
+            total_length = 0
             for x in L:
                 if x.triangle_id % 2 == 0:
@@ -1302 +1303 @@
                 assert x.triangle_id in intersected_triangles
 
+                total_length += x.length
+
+            msg = 'Segments do not add up'
+            assert allclose(total_length, 2), msg
+            
 
     def test_get_intersecting_segments_coinciding(self):
@@ -1333 +1339 @@
 
         # Check all normals point straight down
+        total_length = 0
         for i, x in enumerate(L): 
             assert allclose(x.length, 1.0)
@@ -1345 +1352 @@
             assert x.triangle_id in intersected_triangles
 
-            
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 2), msg            
+
+
 
     def test_get_intersecting_segments2(self):
         """test_get_intersecting_segments(self):
 
-        More complex mesh
-        
+        Lines with a slope
         """
 
-        # Build test mesh (from bounding polygon tests
-        
-
-        
-        pass
+        s2 = sqrt(2.0)/2
+        
+
+        # 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)
+        
+
+        # Diagonal cutting through a vertex and hypothenuses
+        line = [[0, 2], [2, 0]]
+        intersected_triangles = [3,2,5,4]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 4
+
+        #print L
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.length, s2)
+            assert allclose(x.normal, [-s2, -s2])
+            assert allclose(sum(x.normal**2), 1)
+            
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 4*s2), msg
+
+
+        # Diagonal cutting through a vertex and hypothenuses (reversed)
+        line = [[2, 0], [0, 2]]
+        intersected_triangles = [3,2,5,4]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 4
+
+        #print L
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.length, s2)
+            assert allclose(x.normal, [s2, s2])
+            assert allclose(sum(x.normal**2), 1)
+            
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 4*s2), msg                    
+
+
+
+        # Diagonal coinciding with hypothenuses
+        line = [[2, 2], [0, 0]]
+        intersected_triangles = [6,0]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 2
+
+        #print L
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.length, 2*s2)
+            assert allclose(x.normal, [-s2, s2])
+            assert allclose(sum(x.normal**2), 1)
+            
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 4*s2), msg                        
+
+
+        # Diagonal coinciding with hypothenuses (reversed)
+        line = [[0, 0], [2, 2]]
+        intersected_triangles = [6,0]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 2
+
+        #print L
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.length, 2*s2)
+            assert allclose(x.normal, [s2, -s2])
+            assert allclose(sum(x.normal**2), 1)
+            
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 4*s2), msg                        
+
+
+
+        # line with slope [1, -1] cutting through vertices of tri 7 and 6
+        line = [[1, 2], [2, 1]]
+        intersected_triangles = [7,6]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 2
+
+        #print L
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.length, s2)
+            assert allclose(x.normal, [-s2, -s2])
+            assert allclose(sum(x.normal**2), 1)
+            
+            assert x.triangle_id in intersected_triangles
+
+            total_length += x.length
+
+        msg = 'Segments do not add up'
+        assert allclose(total_length, 2*s2), msg
+
+
+        # Arbitrary line with slope [1, -1] cutting through tri 7 and 6
+        line = [[1.1, 2], [2.1, 1]]
+        intersected_triangles = [7,6]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 2
+        
+        # Check all segments
+        total_length = 0
+        for i, x in enumerate(L): 
+            assert allclose(x.normal, [-s2, -s2])
+            assert allclose(sum(x.normal**2), 1)
+
+            msg = 'Triangle %d' %x.triangle_id + ' is not in %s' %(intersected_triangles)
+            assert x.triangle_id in intersected_triangles, msg
+            
+
+
+    def test_get_intersecting_segments3(self):
+        """test_get_intersecting_segments(self):
+
+        Check that line can stop inside a triangle
+        
+        """
+
+
+
+        s2 = sqrt(2.0)/2
+        
+
+        # 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)
+        
+
+        # Line cutting through one triangle and ending on its edge
+        line = [[0.5, 3], [0.5, 1.5]]
+        intersected_triangles = [3]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 1
+        assert L[0].triangle_id == 3
+        assert allclose(L[0].length, 0.5)        
+        assert allclose(L[0].normal, [-1,0])                
+
+
+
+        # Now try to shorten it so that its endpoint falls short of the far edge
+        line = [[0.5, 3], [0.5, 1.6]]
+        intersected_triangles = [3]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 1
+        assert L[0].triangle_id == 3
+        assert allclose(L[0].length, 0.4)
+        assert allclose(L[0].normal, [-1,0])
+
+        intersected_triangles = [3]
+
+        # Now the same, but with direction changed
+        line = [[0.5, 3], [0.5, 1.6]]
+        line = [[0.5, 1.6], [0.5, 3]]        
+        intersected_triangles = [3]                
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 1
+        assert L[0].triangle_id == 3
+        assert allclose(L[0].length, 0.4)
+        assert allclose(L[0].normal, [1,0])                
+        
+
+            
+
+    def test_get_intersecting_segments4(self):
+        """test_get_intersecting_segments(self):
+
+        Simple poly line
+        
+        """
+
+
+
+        s2 = sqrt(2.0)/2
+        
+
+        # 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)
+        
+
+        # Polyline with three segments cutting through mesh
+        line = [[0.5, 3], [0.5, 1.5], [1,1]]
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 2
+
+        for x in L:
+            if x.triangle_id == 3:
+                assert allclose(x.length, 0.5)        
+                assert allclose(x.normal, [-1,0])
+                
+            if x.triangle_id == 2:
+                assert allclose(x.length, s2)
+                assert allclose(x.normal, [-s2,-s2])
+
+
+
+    def test_get_intersecting_segments5(self):
+        """test_get_intersecting_segments(self):
+
+        More complex poly line
+        
+        """
+
+
+
+        s2 = sqrt(2.0)/2
+        
+
+        # 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)
+        
+
+        # Polyline with three segments cutting through mesh
+        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75]] 
+
+        L = mesh.get_intersecting_segments(line)
+        assert len(L) == 3
+
+        for x in L:
+            if x.triangle_id == 3:
+                assert allclose(x.length, 0.5)        
+                assert allclose(x.normal, [-1,0])
+                
+            if x.triangle_id == 2:
+                msg = str(x.length)
+                assert allclose(x.length, s2), msg
+                assert allclose(x.normal, [-s2,-s2])
+
+            if x.triangle_id == 5:
+                segvec = array([line[2][0]-1,
+                                line[2][1]-1])
+                msg = str(x.length)
+                assert allclose(x.length, sqrt(sum(segvec**2))), msg
+                assert allclose(x.normal, [-s2,-s2])                                                
+
+
+    def test_get_intersecting_segments6(self):
+        """test_get_intersecting_segments(self):
+
+        Even more complex poly line, where line breaks within triangle 5
+
+        5 segments are returned even though only four triangles [3,2,5,6] are touched.
+        Triangle 5 therefor has two segments in it.
+        
+        """
+
+
+
+        s2 = sqrt(2.0)/2
+        
+
+        # 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)
+        
+
+        # Polyline with three segments cutting through mesh
+        line = [[0.5, 3], [0.5, 1.5], [1.25, 0.75], [2.25, 1.75]]
+
+        L = mesh.get_intersecting_segments(line)
+        #for x in L:
+        #    print x
+
+        assert len(L) == 5
+
+        for x in L:
+            if x.triangle_id == 3:
+                assert allclose(x.length, 0.5)        
+                assert allclose(x.normal, [-1,0])
+                
+            if x.triangle_id == 2:
+                msg = str(x.length)
+                assert allclose(x.length, s2), msg
+                assert allclose(x.normal, [-s2,-s2])
+
+            if x.triangle_id == 5:
+                if x.segment == ((1.0, 1.0), (1.25, 0.75)):                    
+                    segvec = array([line[2][0]-1,
+                                    line[2][1]-1])
+                    msg = str(x.length)
+                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
+                    assert allclose(x.normal, [-s2,-s2])
+                elif x.segment == ((1.25, 0.75), (1.5, 1.0)):
+                    segvec = array([1.5-line[2][0],
+                                    1.0-line[2][1]])
+                    
+                    assert allclose(x.length, sqrt(sum(segvec**2))), msg
+                    assert allclose(x.normal, [s2,-s2])
+                else:
+                    msg = 'Unknow segment: %s' %x.segment
+                    raise Exception, msg
+                
+
+                    
+            if x.triangle_id == 6:
+                assert allclose(x.normal, [s2,-s2])
+                assert allclose(x.segment, ((1.5, 1.0), (2, 1.5)))
+
+
+      # Internal test that sum of line segments add up
+      # to length of input line
+      #
+      # Could be useful perhaps
+      #
+      #xi1 = line[1][0]
+      #eta1 = line[1][1]
+      #linevector = array([xi1-xi0, eta1-eta0])
+      #linelength = sqrt(sum(linevector**2))
+      #
+      #segmentlength = 0
+      #for segment in triangle_intersections:
+      #    vector = array([segment[1][0] - segment[0][0],
+      #                    segment[1][1] - segment[0][1]])
+      #    length = sqrt(sum(vector**2))      
+      #    segmentlength += length
+      #
+      #msg = 'Sum of intersecting segments do not add up'    
+      #assert allclose(segmentlength, linelength), msg    
+
+
         
 #-------------------------------------------------------------
@@ -1364 +1752 @@
     #suite = unittest.makeSuite(Test_Mesh,'test_two_triangles')
     #suite = unittest.makeSuite(Test_Mesh,'test_get_intersecting_segments_coinciding')
-    suite = unittest.makeSuite(Test_Mesh,'test')
+    suite = unittest.makeSuite(Test_Mesh,'test') #_get_intersecting_segments6')
     runner = unittest.TextTestRunner()
     runner.run(suite)