Changeset 3684

Timestamp:

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

Author:

ole

Message:

Work on get_boundary_polygon using (discontinuous) example that fails

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/general_mesh.py (modified) (7 diffs)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (6 diffs)
• abstract_2d_finite_volumes/test_neighbour_mesh.py (modified) (2 diffs)
• fit_interpolate/interpolate.py (modified) (1 diff)
• utilities/numerical_tools.py (modified) (4 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -80 +80 @@
             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'\
@@ -95 +95 @@
 
 
-        #Register number of elements (N)
+        # Register number of elements (N)
         self.number_of_elements = N = self.triangles.shape[0]
 
@@ -106 +106 @@
 
 
-        #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'
@@ -127 +127 @@
             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
 
@@ -135 +135 @@
 
 
-            #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])
@@ -154 +154 @@
             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()
@@ -178 +178 @@
 
     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

@@ -175 +175 @@
 
     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))
 
@@ -491 +491 @@
 
 
-
-            
         #Start with smallest point and follow boundary (counter clock wise)
         polygon = [p0]      # Storage for final boundary polygon
@@ -508 +506 @@
             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.
@@ -521 +520 @@
                 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
@@ -527 +529 @@
                     # 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    
@@ -536 +536 @@
                 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

@@ -13 +13 @@
 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):
@@ -909 +910 @@
 
 
+    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

@@ -230 +230 @@
         
         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

@@ -4 +4 @@
 """
 
+from math import acos, pi, sqrt
+from warnings import warn
 
 #Establish which Numeric package to use
@@ -21 +23 @@
 
 
+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.
@@ -29 +69 @@
     The angle is measured as a number in [0, 2pi] from v2 to v1.
     """
-    from math import acos, pi, sqrt
   
     # Prepare two Numeric vectors
@@ -41 +80 @@
     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)