Changeset 2709

Ignore:
Timestamp:
Apr 12, 2006, 2:28:05 PM (18 years ago)
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:
6 edited

Unmodified
Removed
• inundation/pyvolution/domain.py

 r2648 """Class Domain - 2D triangular domains for finite-volume computations of the shallow water wave equation conservation laws.
• inundation/pyvolution/mesh.py

 r2683 from general_mesh import General_mesh from math import pi, sqrt class Mesh(General_mesh): #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) 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 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 # 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 from config import epsilon from math import pi from utilities.numerical_tools import anglediff #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

 r2532 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

 r2620 #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 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

 r2704 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

 r2704 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)
Note: See TracChangeset for help on using the changeset viewer.