Changeset 2709 for inundation/pyvolution/mesh.py

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.

File:

• inundation/pyvolution/mesh.py (modified) (7 diffs)

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])