Changeset 2709

Timestamp:

Apr 12, 2006, 2:28:05 PM (19 years ago)

Author:

ole

Message:

Added algorithm for boundary_polygon in the presence of
multiple vertex coordinates. Also tested that new fit_interpolate
can use this new version.

Location:

inundation

Files:

• pyvolution/domain.py (modified) (1 diff)
• pyvolution/mesh.py (modified) (7 diffs)
• pyvolution/test_mesh.py (modified) (1 diff)
• pyvolution/test_shallow_water.py (modified) (2 diffs)
• utilities/numerical_tools.py (modified) (1 diff)
• utilities/test_numerical_tools.py (modified) (1 diff)

inundation/pyvolution/domain.py

@@ -1 +1 @@
 """Class Domain - 2D triangular domains for finite-volume computations of
-   the shallow water wave equation
+   conservation laws.
 
 

inundation/pyvolution/mesh.py

@@ -6 +6 @@
 
 from general_mesh import General_mesh
+from math import pi, sqrt
+        
 
 class Mesh(General_mesh):
@@ -125 +127 @@
                 #of inscribed circle is computed
 
-                from math import sqrt
                 a = sqrt((x0-x1)**2+(y0-y1)**2)
                 b = sqrt((x1-x2)**2+(y1-y2)**2)
@@ -165 +166 @@
     def set_to_inscribed_circle(self,safety_factor = 1):
         #FIXME phase out eventually
-        from math import sqrt
         N = self.number_of_elements
         V = self.vertex_coordinates
@@ -414 +414 @@
 
 
-    def get_boundary_polygon(self):
-        """Return bounding polygon as a list of points
-
-        FIXME: If triangles are listed as discontinuous
-        (e.g vertex coordinates listed multiple times),
-        this may not work as expected.
-        """
+    def get_boundary_polygon(self, verbose = False):
+        """Return bounding polygon for mesh (counter clockwise)
+
+        Using the mesh boundary, derive a bounding polygon for this mesh.
+        If multiple vertex values are present, the algorithm will select the
+        path that contains the mesh.
+        """
+        
         from Numeric import allclose, sqrt, array, minimum, maximum
-
-
-
-        #V = self.get_vertex_coordinates()
+        from utilities.numerical_tools import angle, ensure_numeric        
+
+
+        # Get mesh extent
+        xmin, xmax, ymin, ymax = self.get_extent()
+        pmin = ensure_numeric([xmin, ymin])
+        pmax = ensure_numeric([xmax, ymax])        
+
+
+        # Assemble dictionary of boundary segments and choose starting point 
         segments = {}
-
-        #pmin = (min(self.coordinates[:,0]), min(self.coordinates[:,1]))
-        #pmax = (max(self.coordinates[:,0]), max(self.coordinates[:,1]))
-
-        #FIXME:Can this be written more compactly, e.g.
-        #using minimum and maximium?
-        pmin = array( [min(self.coordinates[:,0]),
-                       min(self.coordinates[:,1]) ] )
-
-        pmax = array( [max(self.coordinates[:,0]),
-                       max(self.coordinates[:,1]) ] )
-
-        mindist = sqrt(sum( (pmax-pmin)**2 ))
+        inverse_segments = {}
+        mindist = sqrt(sum((pmax-pmin)**2)) #Start value across entire mesh
         for i, edge_id in self.boundary.keys():
-            #Find vertex ids for boundary segment
+            # 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))
-
-            #Take the point closest to pmin as starting point
-            #Note: Could be arbitrary, but nice to have
-            #a unique way of selecting
-            dist_A = sqrt(sum( (A-pmin)**2 ))
-            dist_B = sqrt(sum( (B-pmin)**2 ))
-
-            #Find minimal point
+            A = self.get_vertex_coordinate(i, a) # Start
+            B = self.get_vertex_coordinate(i, b) # End
+
+            # Take the point closest to pmin as starting point
+            # Note: Could be arbitrary, but nice to have
+            # a unique way of selecting
+            dist_A = sqrt(sum((A-pmin)**2))
+            dist_B = sqrt(sum((B-pmin)**2))
+
+            #Find lower leftmost point
             if dist_A < mindist:
                 mindist = dist_A
@@ -464 +460 @@
 
 
+            # Sanity check
             if p0 is None:
-                raise 'Weird'
-                p0 = A #We need a starting point (FIXME)
-
-                print 'A', A
-                print 'B', B
-                print 'pmin', pmin
-                print
-
-            segments[A] = B
+                raise Exception('Impossible')
+
+
+            # Register potential paths from A to B
+            if not segments.has_key(tuple(A)):
+                segments[tuple(A)] = [] # Empty list for candidate points                
+                
+            segments[tuple(A)].append(B)                
+
+
 
 
         #Start with smallest point and follow boundary (counter clock wise)
-        polygon = [p0]
-        while len(polygon) < len(self.boundary):
-            p1 = segments[p0]
+        polygon = [p0]      # Storage for final boundary polygon
+        point_registry = {} # Keep track of storage to avoid multiple runs around boundary
+                            # This will only be the case if there are more than one candidate
+                            # FIXME (Ole): Perhaps we can do away with polygon and use
+                            # only point_registry to save space.
+
+        point_registry[tuple(p0)] = 0                            
+                            
+        #while len(polygon) < len(self.boundary):
+        while len(point_registry) < len(self.boundary):
+
+            candidate_list = segments[tuple(p0)]
+            if len(candidate_list) > 1:
+                # Multiple points detected
+                # Take the candidate that is furthest to the clockwise direction,
+                # as that will follow the boundary.
+
+
+                if verbose:
+                    print 'Point %s has multiple candidates: %s' %(str(p0), candidate_list)
+
+
+                # Choose vector against which all angles will be measured
+                if len(polygon) > 1:    
+                    v_prev = p0 - polygon[-2] # Vector that leads to p0
+                else:
+                    # FIXME (Ole): What do we do if the first point has multiple candidates?
+                    # Being the lower left corner, perhaps we can use the
+                    # vector [1, 0], but I really don't know if this is completely
+                    # watertight.
+                    # Another option might be v_prev = [1.0, 0.0]
+                    v_prev = [1.0, 0.0]
+                    
+
+                    
+                # Choose candidate with minimum angle    
+                minimum_angle = 2*pi
+                for pc in candidate_list:
+                    vc = pc-p0  # Candidate vector
+                    
+                    # Angle between each candidate and the previous vector in [-pi, pi]
+                    ac = angle(vc, v_prev)
+                    if ac > pi: ac = pi-ac
+
+                    # take the minimal angle corresponding to the rightmost vector
+                    if ac < minimum_angle:
+                        minimum_angle = ac
+                        p1 = pc             # Best candidate 
+                        
+
+                if verbose is True:
+                    print '  Best candidate %s, angle %f' %(p1, minimum_angle*180/pi)
+                
+            else:
+                p1 = candidate_list[0]
+
+            if point_registry.has_key(tuple(p1)):
+                # We have completed the boundary polygon - yeehaa
+                break 
+            else:
+                point_registry[tuple(p1)] = len(point_registry)
+            
             polygon.append(p1)
             p0 = p1
+
 
         return polygon
@@ -494 +552 @@
 
         from config import epsilon
-        from math import pi
         from utilities.numerical_tools import anglediff
 
@@ -537 +594 @@
 
                 #Normalise
-                from math import sqrt
                 l_u = sqrt(u[0]*u[0] + u[1]*u[1])
                 l_v = sqrt(v[0]*v[0] + v[1]*v[1])                

inundation/pyvolution/test_mesh.py

@@ -849 +849 @@
 
 
+    def test_boundary_polygon_V(self):
+        """Create a discontinuous mesh (duplicate vertices)
+        and check that boundary is as expected
+        
+        """
+        from mesh import Mesh
+        from Numeric import zeros, Float
+        
+
+        #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
+
+        #Duplicate points for triangles edg [4,3,6] (central) and
+        #gid [6,8,7] (top right boundary) to them disconnected
+        #from the others
+
+        e0 = [0.5, 0.5] #9 
+        d0 = [0.5, 0.0] #10
+        g0 = [1.0, 0.5] #11
+        i0 = [1.5, 0.5] #12        
+        
+
+        points = [a, b, c, d, e, f, g, h, i, e0, d0, g0, i0]
+
+
+
+        #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]]
+
+
+        #dea, bae, bec, fgd,
+        #e0d0g0, ghe, gfi, g0i0h
+        vertices = [[3,4,0], [1,0,4], [1,4,2], [5,6,3],
+                    [9,10,11], [6,7,4], [6,5,8], [11,12,7]]        
+
+        mesh = Mesh(points, vertices)
+
+        mesh.check_integrity()
+
+        P = mesh.get_boundary_polygon()
+
+        #print P
+        
+        assert len(P) == 8
+        assert allclose(P, [a, d, f, i, h, e, c, b])
+        assert allclose(P, [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (1.5, 0.5), (1.0, 1.0), (0.5, 0.5), (0.0, 1.0), (0.0, 0.5)])
+        
+
+        for p in points:
+            #print p, P
+            assert inside_polygon(p, P)
 
 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Mesh,'test')
+    suite = unittest.makeSuite(Test_Mesh,'test_boundary_polygon')
+    #suite = unittest.makeSuite(Test_Mesh,'test')
     runner = unittest.TextTestRunner()
     runner.run(suite)

inundation/pyvolution/test_shallow_water.py

@@ -3040 +3040 @@
 
 
-        #Create an triangle shaped domain (reusing coordinates from domain 1),
+        #Create a triangle shaped domain (reusing coordinates from domain 1),
         #formed from the lower and right hand  boundaries and
         #the sw-ne diagonal
@@ -3376 +3376 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Shallow_Water,'test')
-    runner = unittest.TextTestRunner()
+    runner = unittest.TextTestRunner(verbosity=1)    
     runner.run(suite)

inundation/utilities/numerical_tools.py

@@ -57 +57 @@
     
 def anglediff(v0, v1):
-    """Compute difference between angle of vector x0, y0 and x1, y1.
+    """Compute difference between angle of vector v0 (x0, y0) and v1 (x1, y1).
     This is used for determining the ordering of vertices,
     e.g. for checking if they are counter clockwise.

inundation/utilities/test_numerical_tools.py

@@ -47 +47 @@
         
         assert allclose(angle([1.0, -1.0], [1.0, 1.0])/pi*180, 270.0)   
-        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)    
+        assert allclose(angle([1.0, 0.0], [0.0, 1.0])/pi*180, 270.0)
+
+        #From test_get_boundary_polygon_V
+        v_prev = [-0.5, -0.5]
+        vc = [ 0.0,  -0.5]
+        assert allclose(angle(vc, v_prev)/pi*180, 45.0)
+
+        vc = [ 0.5,  0.0]
+        assert allclose(angle(vc, v_prev)/pi*180, 135.0)
+
+        vc = [ -0.5,  0.5]
+        assert allclose(angle(vc, v_prev)/pi*180, 270.0)                
+
+