# Changeset 3684

Ignore:
Timestamp:
Oct 3, 2006, 5:47:08 PM (17 years ago)
Message:

Work on get_boundary_polygon using (discontinuous) example that fails

Location:
anuga_core/source/anuga
Files:
5 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

 r3624 self.geo_reference = geo_reference #Input checks # Input checks msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples. ' msg += 'The supplied array has the shape: %s'\ #Register number of elements (N) # Register number of elements (N) self.number_of_elements = N = self.triangles.shape[0] #Allocate space for geometric quantities # Allocate space for geometric quantities self.normals = zeros((N, 6), Float) self.areas = zeros(N, Float) self.edgelengths = zeros((N, 3), Float) #Get x,y coordinates for all triangles and store # Get x,y coordinates for all triangles and store self.vertex_coordinates = V = self.compute_vertex_coordinates() #Initialise each triangle # Initialise each triangle if verbose: print 'General_mesh: Computing areas, normals and edgelenghts' x2 = V[i, 4]; y2 = V[i, 5] #Area # Area self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2 #Normals #The normal vectors # - point outward from each edge # - are orthogonal to the edge # - have unit length # - Are enumerated according to the opposite corner: #   (First normal is associated with the edge opposite #    the first vertex, etc) # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle # Normals # The normal vectors #   - point outward from each edge #   - are orthogonal to the edge #   - have unit length #   - Are enumerated according to the opposite corner: #     (First normal is associated with the edge opposite #     the first vertex, etc) #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle n0 = array([x2 - x1, y2 - y1]) l2 = sqrt(sum(n2**2)) #Normalise # Normalise n0 /= l0 n1 /= l1 n2 /= l2 #Compute and store # Compute and store self.normals[i, :] = [n0[1], -n0[0], n1[1], -n1[0], n2[1], -n2[0]] #Edgelengths # Edgelengths self.edgelengths[i, :] = [l0, l1, l2] #Build vertex list # Build vertex list if verbose: print 'Building vertex list' self.build_vertexlist() def __repr__(self): return 'Mesh: %d triangles, %d elements'\ return 'Mesh: %d vertex coordinates, %d triangles'\ %(self.coordinates.shape[0], len(self))
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

 r3616 def __repr__(self): return 'Mesh: %d triangles, %d elements, %d boundary segments'\ return 'Mesh: %d vertex_coordinates, %d triangles, %d boundary segments'\ %(self.coordinates.shape[0], len(self), len(self.boundary)) #Start with smallest point and follow boundary (counter clock wise) polygon = [p0]      # Storage for final boundary polygon candidate_list = segments[tuple(p0)] if len(candidate_list) > 1: # Multiple points detected # Multiple points detected (this will be the case for meshes with # duplicate points as those used for discontinuous triangles). # Take the candidate that is furthest to the clockwise direction, # as that will follow the boundary. if len(polygon) > 1: v_prev = p0 - polygon[-2] # Vector that leads to p0 # Normalise #v_prev /= sqrt(sum(v_prev**2)) else: # FIXME (Ole): What do we do if the first point has multiple # 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 for pc in candidate_list: vc = pc-p0  # Candidate vector # Normalise #vc /= sqrt(sum(vc**2)) # Angle between each candidate and the previous vector # in [-pi, pi] ac = angle(vc, v_prev) if ac > pi: ac = pi-ac
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

 r3560 from mesh_factory import rectangular from anuga.config import epsilon from Numeric import allclose, array from Numeric import allclose, array, Int from anuga.coordinate_transforms.geo_reference import Geo_reference from anuga.utilities.polygon import is_inside_polygon from anuga.utilities.numerical_tools import ensure_numeric def distance(x, y): def xtest_boundary_polygon_VI(self): """test_boundary_polygon_VI(self) Create a discontinuous mesh (duplicate vertices) from a real situation that failed and check that boundary is as expected """ from anuga.utilities.polygon import plot_polygons # First do the continuous version of mesh points = [[   6626.85400391,      0.        ], [      0.        ,  38246.4140625 ], [   9656.2734375 ,  68351.265625  ], [  20827.25585938,  77818.203125  ], [  32755.59375   ,  58126.9765625 ], [  35406.3359375 ,  79332.9140625 ], [  31998.23828125,  88799.84375   ], [  23288.65820313, 104704.296875  ], [  32187.57617188, 109816.4375    ], [  50364.08984375, 110763.1328125 ], [  80468.9453125 ,  96184.0546875 ], [  86149.1015625 , 129886.34375   ], [ 118715.359375  , 129886.34375   ], [ 117768.6640625 ,  85770.4296875 ], [ 101485.5390625 ,  45251.9453125 ], [  49985.4140625 ,   2272.06396484], [  51737.94140625,  90559.2109375 ], [  56659.0703125 ,  65907.6796875 ], [  75735.4765625 ,  23762.00585938], [  52341.70703125,  38563.39453125]] triangles = [[19, 0,15], [ 2, 4, 3], [ 4, 2, 1], [ 1,19, 4], [15,18,19], [18,14,17], [19, 1, 0], [ 6, 8, 7], [ 8, 6,16], [10, 9,16], [17, 5, 4], [16,17,10], [17,19,18], [ 5,17,16], [10,14,13], [10,17,14], [ 8,16, 9], [12,11,10], [10,13,12], [19,17, 4], [16, 6, 5]] mesh = Mesh(points, triangles) mesh.check_integrity() Pref = mesh.get_boundary_polygon() plot_polygons([ensure_numeric(Pref)], 'goodP') for p in points: assert is_inside_polygon(p, Pref) # Then do the discontinuous version import warnings warnings.filterwarnings('ignore') points = [[  52341.70703125,  38563.39453125], [   6626.85400391,      0.        ], [  49985.4140625 ,   2272.06396484], [   9656.2734375 ,  68351.265625  ], [  32755.59375   ,  58126.9765625 ], [  20827.25585938,  77818.203125  ], [  32755.59375   ,  58126.9765625 ], [   9656.2734375 ,  68351.265625  ], [      0.        ,  38246.4140625 ], [      0.        ,  38246.4140625 ], [  52341.70703125,  38563.39453125], [  32755.59375   ,  58126.9765625 ], [  49985.4140625 ,   2272.06396484], [  75735.4765625 ,  23762.00585938], [  52341.70703125,  38563.39453125], [  75735.4765625 ,  23762.00585938], [ 101485.5390625 ,  45251.9453125 ], [  56659.0703125 ,  65907.6796875 ], [  52341.70703125,  38563.39453125], [      0.        ,  38246.4140625 ], [   6626.85400391,      0.        ], [  31998.23828125,  88799.84375   ], [  32187.57617188, 109816.4375    ], [  23288.65820313, 104704.296875  ], [  32187.57617188, 109816.4375    ], [  31998.23828125,  88799.84375   ], [  51737.94140625,  90559.2109375 ], [  80468.9453125 ,  96184.0546875 ], [  50364.08984375, 110763.1328125 ], [  51737.94140625,  90559.2109375 ], [  56659.0703125 ,  65907.6796875 ], [  35406.3359375 ,  79332.9140625 ], [  32755.59375   ,  58126.9765625 ], [  51737.94140625,  90559.2109375 ], [  56659.0703125 ,  65907.6796875 ], [  80468.9453125 ,  96184.0546875 ], [  56659.0703125 ,  65907.6796875 ], [  52341.70703125,  38563.39453125], [  75735.4765625 ,  23762.00585938], [  35406.3359375 ,  79332.9140625 ], [  56659.0703125 ,  65907.6796875 ], [  51737.94140625,  90559.2109375 ], [  80468.9453125 ,  96184.0546875 ], [ 101485.5390625 ,  45251.9453125 ], [ 117768.6640625 ,  85770.4296875 ], [  80468.9453125 ,  96184.0546875 ], [  56659.0703125 ,  65907.6796875 ], [ 101485.5390625 ,  45251.9453125 ], [  32187.57617188, 109816.4375    ], [  51737.94140625,  90559.2109375 ], [  50364.08984375, 110763.1328125 ], [ 118715.359375  , 129886.34375   ], [  86149.1015625 , 129886.34375   ], [  80468.9453125 ,  96184.0546875 ], [  80468.9453125 ,  96184.0546875 ], [ 117768.6640625 ,  85770.4296875 ], [ 118715.359375  , 129886.34375   ], [  52341.70703125,  38563.39453125], [  56659.0703125 ,  65907.6796875 ], [  32755.59375   ,  58126.9765625 ], [  51737.94140625,  90559.2109375 ], [  31998.23828125,  88799.84375   ], [  35406.3359375 ,  79332.9140625 ]] #points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation triangles = [[ 0, 1, 2], [ 3, 4, 5], [ 6, 7, 8], [ 9,10,11], [12,13,14], [15,16,17], [18,19,20], [21,22,23], [24,25,26], [27,28,29], [30,31,32], [33,34,35], [36,37,38], [39,40,41], [42,43,44], [45,46,47], [48,49,50], [51,52,53], [54,55,56], [57,58,59], [60,61,62]] mesh = Mesh(points, triangles) mesh.check_integrity() P = mesh.get_boundary_polygon() plot_polygons([ensure_numeric(P)], 'badP') print P for p in points: #assert is_inside_polygon(p, P) if not is_inside_polygon(p, P): print 'Point %s is not in P' %(p) else: print 'Point %s OK' %(p) #for i in range(len(Pref)): #    print P[i], Pref[i] # #print ensure_numeric(P) - ensure_numeric(Pref) def test_lone_vertices(self): a = [2.0, 1.0]
• ## anuga_core/source/anuga/fit_interpolate/interpolate.py

 r3633 if verbose: print 'Getting indices inside mesh boundary' #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() self.inside_poly_indices, self.outside_poly_indices  = \ in_and_outside_polygon(point_coordinates, self.mesh.get_boundary_polygon(), closed = True, verbose = verbose) in_and_outside_polygon(point_coordinates, self.mesh.get_boundary_polygon(), closed = True, verbose = verbose) #print "self.inside_poly_indices",self.inside_poly_indices #print "self.outside_poly_indices",self.outside_poly_indices
• ## anuga_core/source/anuga/utilities/numerical_tools.py

 r3676 """ from math import acos, pi, sqrt from warnings import warn #Establish which Numeric package to use def safe_acos(x): """Safely compute acos Protect against cases where input argument x is outside the allowed interval [-1.0, 1.0] by no more than machine precision """ error_msg = 'Input to acos is outside allowed domain [-1.0, 1.0]. I got %.12f' %x warning_msg = 'Changing argument to acos from %.18f to %.1f' %(x, sign(x)) # FIXME(Ole): Need function to compute machine precision as this is also # used in fit_interpolate/search_functions.py eps = 1.0e-15  # Machine precision if x < -1.0: if x < -1.0 - eps: raise ValueError, errmsg else: warn(warning_msg) x = -1.0 if x > 1.0: if x > 1.0 + eps: raise ValueError, errmsg else: print 'NOTE: changing argument to acos from %.18f to 1.0' %x x = 1.0 return acos(x) def sign(x): if x > 0: return 1 if x < 0: return -1 if x == 0: return 0 def angle(v1, v2=None): """Compute angle between 2D vectors v1 and v2. The angle is measured as a number in [0, 2pi] from v2 to v1. """ from math import acos, pi, sqrt # Prepare two Numeric vectors v1 = v1/sqrt(sum(v1**2)) v2 = v2/sqrt(sum(v2**2)) # Compute angle p = innerproduct(v1, v2) c = innerproduct(v1, normal_vector(v2)) # Projection onto normal # (negative cross product) #print "p",p #print "v1", v1 #print "v2", v2 # Warning, this is a hack.  It could cause code to go in loop forever if False: try: theta = acos(p) #print "theta",theta except ValueError: print "Doing a hack in numerical tools." print "p",p print "v1", v1 print "v2", v2 if p > (1.0 - 1e-12): #sus, checking a float # Throw a warning theta = 0.0 else: raise else: theta = acos(p) theta = safe_acos(p) #   print "problem with p",p # as p goes to 1 theta goes to 0 # Correct if v1 is in quadrant 3 or 4 with respect to v2 (as the x-axis)
Note: See TracChangeset for help on using the changeset viewer.