Changeset 984

Timestamp:

Mar 1, 2005, 6:26:00 PM (20 years ago)

Author:

ole

Message:

Added get_boundary_polygon to mesh and started on a pre-crop of data points within least squares

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (1 diff)
• mesh.py (modified) (2 diffs)
• test_mesh.py (modified) (1 diff)
• test_util.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -395 +395 @@
 
 
+        #Remove points falling outside mesh boundary
+        #from util import inside_polygon
+        #P = self.mesh.get_boundary_polygon()
+        #FIXME: TODO
+
+
 
         

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -48 +48 @@
     Mesh takes the optional third argument boundary which is a
     dictionary mapping from (element_id, edge_id) to boundary tag.
-    The default value is None which will assign the defualt_boundary_tag
+    The default value is None which will assign the default_boundary_tag
     as specified in config.py to all boundary edges.
     """
@@ -346 +346 @@
         return tags.keys()    
             
-    
+
+    def get_boundary_polygon(self):
+        """Return bounding polygon as a list of points
+        """
+        from Numeric import allclose
+        
+
+        V = self.get_vertex_coordinates()
+        segments = {}
+
+        pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
+        for i, edge_id in self.boundary.keys():
+            #Find vertex ids for boundary segment
+            if edge_id == 0: a = 1; b = 2
+            if edge_id == 1: a = 2; b = 0
+            if edge_id == 2: a = 0; b = 1                
+
+            A = tuple(self.get_vertex_coordinate(i, a))
+            B = tuple(self.get_vertex_coordinate(i, b))
+
+            #Find minimal point
+            if allclose(A, pmin): p0 = A
+            if allclose(B, pmin): p0 = B                
+
+            segments[A] = B
+
+        
+        #Start with smallest point and follow boundary (counter clock wise)
+        polygon = [p0]
+        while len(polygon) < len(self.boundary):
+            p1 = segments[p0]
+            polygon.append(p1)
+            p0 = p1
+            
+        return polygon
+        
+            
+
+        
 
     def check_integrity(self):

inundation/ga/storm_surge/pyvolution/test_mesh.py

@@ -649 +649 @@
 
         
-        
+    def test_boundary_polygon(self):
+        from mesh_factory import rectangular
+        from mesh import Mesh
+        from Numeric import zeros, Float
+        from util import inside_polygon
+        
+        #Create basic mesh
+        points, vertices, boundary = rectangular(2, 2)
+        mesh = Mesh(points, vertices, boundary)
+        
+
+        P = mesh.get_boundary_polygon()
+
+        assert len(P) == 8
+        assert allclose(P, [[0.0, 0.0], [0.5, 0.0], [1.0, 0.0],
+                            [1.0, 0.5], [1.0, 1.0], [0.5, 1.0],
+                            [0.0, 1.0], [0.0, 0.5]])
+        for p in points:
+            #print p, P
+            assert inside_polygon(p, P) 
+
+
+    def test_boundary_polygon_II(self):
+        from mesh import Mesh
+        from Numeric import zeros, Float
+        from util import inside_polygon
+        
+        #Points
+        a = [0.0, 0.0] #0
+        b = [0.0, 0.5] #1
+        c = [0.0, 1.0] #2
+        d = [0.5, 0.0] #3
+        e = [0.5, 0.5] #4
+        f = [1.0, 0.0] #5
+        g = [1.0, 0.5] #6
+        h = [1.0, 1.0] #7
+        i = [1.5, 0.5] #8
+
+        points = [a, b, c, d, e, f, g, h, i]
+
+        #dea, bae, bec, fgd,
+        #edg, ghe, gfi, gih
+        vertices = [ [3,4,0], [1,0,4], [1,4,2], [5,6,3],
+                     [4,3,6], [6,7,4], [6,5,8], [6,8,7]]
+
+        mesh = Mesh(points, vertices)
+        
+        mesh.check_integrity()
+        
+        P = mesh.get_boundary_polygon()
+
+        assert len(P) == 8
+        assert allclose(P, [a, d, f, i, h, e, c, b])
+
+        for p in points:
+            #print p, P
+            assert inside_polygon(p, P) 
+
+
+    def test_boundary_polygon_III(self):
+        """Same as II but vertices ordered differently
+        """
+        
+        from mesh import Mesh
+        from Numeric import zeros, Float
+        from util import inside_polygon
+        
+        #Points
+        a = [0.0, 0.0] #0
+        b = [0.0, 0.5] #1
+        c = [0.0, 1.0] #2
+        d = [0.5, 0.0] #3
+        e = [0.5, 0.5] #4
+        f = [1.0, 0.0] #5
+        g = [1.0, 0.5] #6
+        h = [1.0, 1.0] #7
+        i = [1.5, 0.5] #8
+
+        points = [a, b, c, d, e, f, g, h, i]
+
+        #edg, ghe, gfi, gih
+        #dea, bae, bec, fgd,
+        vertices = [[4,3,6], [6,7,4], [6,5,8], [6,8,7],
+                    [3,4,0], [1,0,4], [1,4,2], [5,6,3]]
+
+
+        mesh = Mesh(points, vertices)
+        mesh.check_integrity()
+
+        
+        P = mesh.get_boundary_polygon()
+
+        assert len(P) == 8
+        assert allclose(P, [a, d, f, i, h, e, c, b])
+
+        for p in points:
+            assert inside_polygon(p, P) 
+
+
         
 #-------------------------------------------------------------

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -866 +866 @@
         assert not inside_polygon( (0.5, -0.5), polygon )               
         
-        #Now try the vector formulation returning indices
-        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        res = inside_polygon( points, polygon )
-        
-        assert allclose( res, [0,1,2] )
         
         #Very convoluted polygon
@@ -913 +908 @@
         assert inside_polygon( (0, -0.5), polygon )
         assert not inside_polygon( (0.5, 0.5), polygon )        
-        
+
+    def test_inside_polygon_vector_version(self):       
+        #Now try the vector formulation returning indices
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]        
+        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        res = inside_polygon( points, polygon )
+        
+        assert allclose( res, [0,1,2] )
+
 
     def test_populate_polygon(self):