Changeset 9242

Timestamp:

Jul 4, 2014, 11:31:31 AM (10 years ago)

Author:

davies

Message:

Improvements + tests for polygon intersections in spatialInputUtils

Location:

trunk/anuga_core/source/anuga/utilities

Files:

• spatialInputUtil.py (modified) (10 diffs)
• test_spatialInputUtil.py (modified) (1 diff)

trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py

@@ -352 +352 @@
     
     #################################################################################
-    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold):
+    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold, verbose=False):
         """
             Add intersectionPt to line_pts, either by inserting it, or if a point on line_pts is
@@ -381 +381 @@
             # account for the possibility of L1_pts having > 2 coordinates
             # (riverWalls)
+            if verbose:
+                print '      Inserting new point'
             dummyPt=copy.copy(L1_pts[tmp[1]])
             L1_pts.insert(tmp[1]+1,dummyPt) 
@@ -390 +392 @@
                 # Set 3rd coordinate as distance-weighted average of the others
                 d0=((L1_pts[tmp[1]][0]-L1_pts[tmp[1]+1][0])**2.+\
-                   (L1_pts[tmp[1]][0]-L1_pts[tmp[1]+1][1])**2.)**0.5
+                   (L1_pts[tmp[1]][1]-L1_pts[tmp[1]+1][1])**2.)**0.5
                 d1=((L1_pts[tmp[1]+2][0]-L1_pts[tmp[1]+1][0])**2.+\
                    (L1_pts[tmp[1]+2][1]-L1_pts[tmp[1]+1][1])**2.)**0.5
@@ -396 +398 @@
     
         else:
+            if verbose:
+                print '      Shifting existing point'
             # Move a point already on L1
             L1_pts=shift_point_on_line(iP, L1_pts, tmp[1])
@@ -435 +439 @@
     #######################################################################################################
     
-    def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0, buf=1.0e-06, tol2 = 100,
+    def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0, buf=1.0e-05, tol2 = 100,
                                   verbose=True, nameFlag=''):
         """
@@ -495 +499 @@
             # Insert the points into the line segments
             for i in range(len(intersectionPts)):
-                L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, point_movement_threshold)
-                L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, point_movement_threshold)
+                L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, point_movement_threshold, verbose=verbose)
+                L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, point_movement_threshold, verbose=verbose)
                 
             # Convert to the input format
@@ -796 +800 @@
                                 verbose=verbose, nameFlag='Bounding Pol intersects '+ n1)
                 riverWalls[n1]=Wkb2ListPts(rw1)
-                bounding_polygon=Wkb2ListPts(bp2,removeLast=True)
+                # Since the bounding polygon is a loop, the first/last points are the same
+                # If one of these was moved, the other should be moved too. Since we
+                # will drop the last bounding_polygon point, we only need to worry about the first
+                bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
+                if(bounding_polygon[-1] is not bounding_polygon[0]):
+                    bounding_polygon[0]=bounding_polygon[-1]
+                # Drop the last point
+                bounding_polygon=bounding_polygon[:-1]
     
         # Clean intersections of bounding polygon and breaklines
@@ -813 +824 @@
                                 verbose=verbose, nameFlag='Bounding Pol intersects '+n1)
                 breakLines[n1]=Wkb2ListPts(bl1)
-                bounding_polygon=Wkb2ListPts(bp2, removeLast=True)
-      
+                # Since the bounding polygon is a loop, the first/last points are the same
+                # If one of these was moved, the other should be moved too. Since we
+                # will drop the last bp2 point, we only need to worry about the first
+                bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
+                if(bounding_polygon[-1] is not bounding_polygon[0]):
+                    bounding_polygon[0]=bounding_polygon[-1]
+                # Drop the last point
+                bounding_polygon=bounding_polygon[:-1]
+
         # Remove the extra 0.0 from bounding polygon [this cannot have 3 coordinates] 
         bounding_polygon = [ bounding_polygon[i][0:2] for i in range(len(bounding_polygon))]
@@ -820 +838 @@
         for n1 in breakLines.keys():
             breakLines[n1] = [breakLines[n1][i][0:2] for i in range(len(breakLines[n1]))]
+
+        # Check that all mesh-boundary points are inside the bounding polygon
+        from anuga.geometry.polygon import outside_polygon
+        for blCat in [riverWalls, breakLines]:
+            for n1 in blCat.keys():
+                l=len(blCat[n1])
+                # Test every point -- means we can strip 3rd coordinate if needed
+                for j in range(l):
+                    isOut=outside_polygon(blCat[n1][j][0:2], bounding_polygon)
+                    if(len(isOut)>0):
+                        msg='Breakline/riverwall point '+str(blCat[n1][j][0:2]) +' on '+ n1+\
+                            ' is outside the bounding polygon.\n'+\
+                            'Check that it exceeds the bounding polygon by a distance < point_movement_threshold \n'+\
+                            ' so it can be moved back onto the polygon'
+                        print 'Polygon\n '
+                        print bounding_polygon
+                        print 'Line \n'
+                        print blCat[n1]
+                        raise Exception, msg
      
         return [bounding_polygon, breakLines, riverWalls]
@@ -961 +998 @@
         raise ImportError, msg
     
-    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold):
+    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold,verbose=False):
         raise ImportError, msg
 

trunk/anuga_core/source/anuga/utilities/test_spatialInputUtil.py

@@ -362 +362 @@
             assert True
 
-        #print 'After first add_iinersections'
         #################################################################
         # Fix the riverwall, and it should work
-        riverWalls={ 'rw1': [[-0.01, 8., 2.],[10.01, 4., 3.]],
-                     'rw2': [[5.3, -0.01, 1.], [10., 10.01, 2.]] }
+        riverWalls={ 'rw1': [[-0.000001, 8., 2.],[10.0000001, 4., 3.]]
+                      }
         # This should work
         newBP, newBL, newRW=su.add_intersections_to_domain_features(bounding_polygon,\
                                 breakLines, riverWalls, point_movement_threshold=0.02,\
                                 verbose=verbose)
+
+        # There should be several new points on breakLines + riverWalls
+        assert newBL['bl1'][1]==newBL['bl2'][1]
+        assert newRW['rw1'][2][0:2]==newBL['bl1'][2]
+        assert newRW['rw1'][1][0:2]==newBL['bl2'][2]
+        # rw1 x/y coords are close to y=8-3/10*x
+        #     x/z coords are close to z=2+x/10
+        assert numpy.allclose(newRW['rw1'][1][2], newRW['rw1'][1][0]/10.+2)
+        assert numpy.allclose(newRW['rw1'][2][2], newRW['rw1'][2][0]/10.+2)
 
         return