Changeset 316

Timestamp:

Sep 16, 2004, 5:31:08 PM (21 years ago)

Author:

duncan

Message:

bigger least squares tests. Making it more stand-alone

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (7 diffs)
• test_least_squares.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -7 +7 @@
 
 
-from mesh import Mesh
-
+"""Classes implementing general 2D geometrical mesh of triangles.
+
+   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Geoscience Australia, 2004   
+"""
+
+
+class Mesh:
+    """Collection of triangular elements (purely geometric)
+
+    A triangular element is defined in terms of three vertex ids,
+    ordered counter clock-wise, 
+    each corresponding to a given coordinate set.
+    Vertices from different elements can point to the same
+    coordinate set.
+
+    Coordinate sets are implemented as an N x 2 Numeric array containing
+    x and y coordinates. 
+    
+
+    To instantiate:
+       Mesh(coordinates, triangles)
+
+    where
+
+      coordinates is either a list of 2-tuples or an Mx2 Numeric array of
+      floats representing all x, y coordinates in the mesh.
+
+      triangles is either a list of 3-tuples or an Nx3 Numeric array of
+      integers representing indices of all vertices in the mesh.
+      Each vertex is identified by its index i in [0, M-1].
+
+
+    Example:
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0,0.0]
+        e = [2.0, 2.0]
+        
+        points = [a, b, c, e]
+        triangles = [ [1,0,2], [1,2,3] ]   #bac, bce 
+        mesh = Mesh(points, triangles)
+    
+        #creates two triangles: bac and bce
+
+        This is a cut down version of mesh from pyvolution\mesh.py
+    """
+
+    def __init__(self, coordinates, triangles):
+        """
+        Build triangles from x,y coordinates (sequence of 2-tuples or
+        Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
+        or Nx3 Numeric array of non-negative integers).
+        """
+
+        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
+
+        self.triangles = array(triangles).astype(Int)
+        self.coordinates = array(coordinates)
+
+        #Input checks
+        msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.triangles.shape) == 2, msg
+
+        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
+        assert len(self.coordinates.shape) == 2, msg        
+
+        msg = 'Vertex indices reference non-existing coordinate sets'
+        assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
+
+
+        #Register number of elements (N)
+        self.number_of_elements = N = self.triangles.shape[0]
+    
+        
+        self.normals = zeros((N, 6), Float)
+
+        #Get x,y coordinates for all triangles and store
+        self.vertex_coordinates = V = self.compute_vertex_coordinates()
+        
+        #Initialise each triangle
+        for i in range(N):
+            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
+            
+            x0 = V[i, 0]; y0 = V[i, 1]
+            x1 = V[i, 2]; y1 = V[i, 3]
+            x2 = V[i, 4]; y2 = V[i, 5]            
+
+            #Area
+            area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
+
+            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
+            msg += ' is degenerate:  area == %f' %area
+            assert area > 0.0, msg
+        
+
+            #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])
+            l0 = sqrt(sum(n0**2))
+
+            n1 = array([x0 - x2, y0 - y2])
+            l1 = sqrt(sum(n1**2))
+
+            n2 = array([x1 - x0, y1 - y0])
+            l2 = sqrt(sum(n2**2))
+
+            #Normalise
+            n0 /= l0
+            n1 /= l1
+            n2 /= l2
+
+            #Compute and store
+            self.normals[i, :] = [n0[1], -n0[0],
+                                  n1[1], -n1[0],
+                                  n2[1], -n2[0]]             
+
+         
+        #Build vertex list
+        self.build_vertexlist()        
+
+
+        
+    def __len__(self):
+        return self.number_of_elements
+
+    def __repr__(self):
+        return 'Mesh: %d triangles, %d elements'\
+               %(self.coordinates.shape[0], len(self))
+
+
+
+    #DSG I don't know if I need this
+    def build_vertexlist(self):
+        pass
+    
+    # maybe add this latter..    
+    def check_integrity(self):
+        pass
+    
+
+    def get_normals(self):
+        """Return all normal vectors.
+        Return normal vectors for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.normals
+
+
+    def get_normal(self, i, j):
+        """Return normal vector j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.normals[i, 2*j:2*j+2]     
+    
+
+            
+    def get_vertex_coordinates(self):
+        """Return all vertex coordinates.
+        Return all vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
+        """
+        return self.vertex_coordinates
+
+
+    def get_vertex_coordinate(self, i, j):
+        """Return coordinates for vertex j of the i'th triangle.
+        Return value is the numeric array slice [x, y]
+        """
+        return self.vertex_coordinates[i, 2*j:2*j+2]     
+
+
+    def compute_vertex_coordinates(self):
+        """Return vertex coordinates for all triangles as an Nx6 array
+        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
+        """
+
+        #FIXME: Perhaps they should be ordered as in obj files??
+        #See quantity.get_vertex_values
+        
+        from Numeric import zeros, Float
+
+        N = self.number_of_elements
+        vertex_coordinates = zeros((N, 6), Float) 
+
+        for i in range(N):
+            for j in range(3):
+                k = self.triangles[i,j]  #Index of vertex 0
+                v_k = self.coordinates[k]
+                vertex_coordinates[i, 2*j+0] = v_k[0]
+                vertex_coordinates[i, 2*j+1] = v_k[1]
+
+        return vertex_coordinates
+  
+    def get_triangles(self, unique=False):
+        """Get connectivity
+        If unique is True give them only once as stored internally.
+        If unique is False return
+        """
+
+        if unique is True:
+            return self.triangles
+        else:
+            from Numeric import reshape, array, Int
+                
+            m = len(self)  #Number of volumes
+            M = 3*m        #Total number of unique vertices
+            return reshape(array(range(M)).astype(Int), (m,3))
+
+
+
+from Numeric import zeros, array, Float, Int, dot,transpose
+from LinearAlgebra import solve_linear_equations
+
+EPSILON = 1.0e-13
 
 class Interpolation(Mesh):
 
-    def __init__(self, vertex_coordinates, triangles, data_coordinates):
+    def __init__(self, vertex_coordinates, triangles, point_coordinates):
         """ Build interpolation matrix mapping from
         function values at vertices to function values at data points
@@ -24 +245 @@
           integers representing indices of all vertices in the mesh.
 
-          data_coordinates: List of coordinate pairs [x, y] of data points
+          point_coordinates: List of coordinate pairs [x, y] of data points
           (or an nx2 Numeric array)
           
         """
 
-        from Numeric import zeros, array, Float, Int, dot
-        from config import epsilon
-
         #Convert input to Numeric arrays
-        data_coordinates = array(data_coordinates).astype(Float)
+        point_coordinates = array(point_coordinates).astype(Float)
         vertex_coordinates = array(vertex_coordinates).astype(Float)
         triangles = array(triangles).astype(Int)                
@@ -43 +261 @@
         #Build n x m interpolation matrix        
         m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex)
-        n = data_coordinates.shape[0]               #Number of data points         
+        n = point_coordinates.shape[0]               #Number of data points         
 
         self.A_m = zeros((n,m), Float)
@@ -51 +269 @@
             #For each data_coordinate point
 
-            x = data_coordinates[i]
+            x = point_coordinates[i]
             for k in range(len(self)):
                 #For each triangle (brute force)
@@ -72 +290 @@
 
                 #FIXME: Maybe move out to test or something
-                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < epsilon
+                assert abs(sigma0 + sigma1 + sigma2 - 1.0) < EPSILON
 
                 
@@ -89 +307 @@
 
     def smooth_to_mesh(self, z):
-        from Numeric import zeros, array, Float,transpose,dot
-        from LinearAlgebra import solve_linear_equations
+        """ Linearly interpolate point data values to the vertices.
+        
+        Inputs:
+          z: Vector of data at the point_coordinates.  One value per point.
+          
+        """
         #Convert input to Numeric arrays
         z = array(z).astype(Float)
@@ -100 +322 @@
         f = solve_linear_equations(self.AtA_m,Atz_m)
         return f
+
+    def smooth_to_points(self, f):
+        """ Linearly interpolate point data values to the vertices.
+        
+        Inputs:
+          z: Vector of data at the point_coordinates.  One value per point.
+          
+        """
+        return dot(self.A_m,f)

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -159 +159 @@
         assert allclose(f, answer)
         
+    def test_smooth_to_meshII(self):
+        a = [0.0, 0.0]
+        b = [0.0, 5.0]
+        c = [5.0, 0.0]
+        points = [a, b, c]
+        triangles = [ [1,0,2] ]   #bac
+         
+        d1 = [1.0, 1.0]
+        d2 = [1.0, 2.0]
+        d3 = [3.0,1.0]
+        data_coords = [ d1, d2, d3]       
+        z = self.linear_function(data_coords)
+        mesh = Interpolation(points, triangles, data_coords)
+        f =  mesh.smooth_to_mesh(z)
+        answer = self.linear_function(points)
+        assert allclose(f, answer)
+    
+    def test_smooth_to_meshIII(self):
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
+
+         
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0,2.0],
+                        [1.0, 1.0],
+                        [2.0, 1.0],
+                        [0.0,0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2,-0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0,1.0]]
+        
+        z = self.linear_function(point_coords)
+        mesh = Interpolation(vertices, triangles, point_coords)
+        f =  mesh.smooth_to_mesh(z)
+        answer = self.linear_function(vertices)
+        assert allclose(f, answer)
+        
+    def linear_function(self,point):
+        point = array(point)
+        return point[:,0]+point[:,1]
+    
+    def test_smooth_to_points(self):
+        v0 = [0.0, 0.0]
+        v1 = [0.0, 5.0]
+        v2 = [5.0, 0.0]
+        
+        vertices = [v0, v1, v2]
+        triangles = [ [1,0,2] ]   #bac
+        
+         
+        d0 = [1.0, 1.0]
+        d1 = [1.0, 2.0]
+        d2 = [3.0,1.0]
+        point_coords = [ d0, d1, d2]
+        
+        mesh = Interpolation(vertices, triangles, point_coords)
+        f = self.linear_function(vertices)
+        z =  mesh.smooth_to_points(f)
+        answer = self.linear_function(point_coords)
+        #print "z",z
+        #print "answer",answer
+        assert allclose(z, answer)
+       
+    def test_smooth_to_pointsII(self):
+        a = [-1.0, 0.0]
+        b = [3.0, 4.0]
+        c = [4.0,1.0]
+        d = [-3.0, 2.0] #3
+        e = [-1.0,-2.0]
+        f = [1.0, -2.0] #5
+
+        vertices = [a, b, c, d,e,f]
+        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
+
+         
+        point_coords = [[-2.0, 2.0],
+                        [-1.0, 1.0],
+                        [0.0,2.0],
+                        [1.0, 1.0],
+                        [2.0, 1.0],
+                        [0.0,0.0],
+                        [1.0, 0.0],
+                        [0.0, -1.0],
+                        [-0.2,-0.5],
+                        [-0.9, -1.5],
+                        [0.5, -1.9],
+                        [3.0,1.0]]
+        
+        mesh = Interpolation(vertices, triangles, point_coords)
+        f = self.linear_function(vertices)
+        z =  mesh.smooth_to_points(f)
+        answer = self.linear_function(point_coords)
+        #print "z",z
+        #print "answer",answer
+        assert allclose(z, answer) 
 #-------------------------------------------------------------
 if __name__ == "__main__":