Changeset 321


Ignore:
Timestamp:
Sep 20, 2004, 9:19:01 AM (21 years ago)
Author:
duncan
Message:

adding a general mesh file, that pyvolution mesh inherits from

Location:
inundation/ga/storm_surge/pyvolution
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r318 r321  
    1313"""
    1414
    15 #FIXME: What is the reason for duplication class Mesh
    16 #(or large parts of it at least) here?
    17 
    18 class Mesh:
    19     """Collection of triangular elements (purely geometric)
    20 
    21     A triangular element is defined in terms of three vertex ids,
    22     ordered counter clock-wise,
    23     each corresponding to a given coordinate set.
    24     Vertices from different elements can point to the same
    25     coordinate set.
    26 
    27     Coordinate sets are implemented as an N x 2 Numeric array containing
    28     x and y coordinates.
    29    
    30 
    31     To instantiate:
    32        Mesh(coordinates, triangles)
    33 
    34     where
    35 
    36       coordinates is either a list of 2-tuples or an Mx2 Numeric array of
    37       floats representing all x, y coordinates in the mesh.
    38 
    39       triangles is either a list of 3-tuples or an Nx3 Numeric array of
    40       integers representing indices of all vertices in the mesh.
    41       Each vertex is identified by its index i in [0, M-1].
     15EPSILON = 1.0e-13
    4216
    4317
    44     Example:
    45         a = [0.0, 0.0]
    46         b = [0.0, 2.0]
    47         c = [2.0,0.0]
    48         e = [2.0, 2.0]
    49        
    50         points = [a, b, c, e]
    51         triangles = [ [1,0,2], [1,2,3] ]   #bac, bce
    52         mesh = Mesh(points, triangles)
    53    
    54         #creates two triangles: bac and bce
    55 
    56         This is a cut down version of mesh from pyvolution\mesh.py
    57     """
    58 
    59     def __init__(self, coordinates, triangles):
    60         """
    61         Build triangles from x,y coordinates (sequence of 2-tuples or
    62         Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
    63         or Nx3 Numeric array of non-negative integers).
    64         """
    65 
    66         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
    67 
    68         self.triangles = array(triangles).astype(Int)
    69         self.coordinates = array(coordinates)
    70 
    71         #Input checks
    72         msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples'
    73         assert len(self.triangles.shape) == 2, msg
    74 
    75         msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
    76         assert len(self.coordinates.shape) == 2, msg       
    77 
    78         msg = 'Vertex indices reference non-existing coordinate sets'
    79         assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
    80 
    81 
    82         #Register number of elements (N)
    83         self.number_of_elements = N = self.triangles.shape[0]
    84    
    85        
    86         self.normals = zeros((N, 6), Float)
    87 
    88         #Get x,y coordinates for all triangles and store
    89         self.vertex_coordinates = V = self.compute_vertex_coordinates()
    90        
    91         #Initialise each triangle
    92         for i in range(N):
    93             #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
    94            
    95             x0 = V[i, 0]; y0 = V[i, 1]
    96             x1 = V[i, 2]; y1 = V[i, 3]
    97             x2 = V[i, 4]; y2 = V[i, 5]           
    98 
    99             #Area
    100             area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
    101 
    102             msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
    103             msg += ' is degenerate:  area == %f' %area
    104             assert area > 0.0, msg
    105        
    106 
    107             #Normals
    108             #The normal vectors
    109             # - point outward from each edge
    110             # - are orthogonal to the edge
    111             # - have unit length
    112             # - Are enumerated according to the opposite corner:
    113             #   (First normal is associated with the edge opposite
    114             #    the first vertex, etc)
    115             # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
    116 
    117             n0 = array([x2 - x1, y2 - y1])
    118             l0 = sqrt(sum(n0**2))
    119 
    120             n1 = array([x0 - x2, y0 - y2])
    121             l1 = sqrt(sum(n1**2))
    122 
    123             n2 = array([x1 - x0, y1 - y0])
    124             l2 = sqrt(sum(n2**2))
    125 
    126             #Normalise
    127             n0 /= l0
    128             n1 /= l1
    129             n2 /= l2
    130 
    131             #Compute and store
    132             self.normals[i, :] = [n0[1], -n0[0],
    133                                   n1[1], -n1[0],
    134                                   n2[1], -n2[0]]             
    135 
    136          
    137         #Build vertex list
    138         self.build_vertexlist()       
    139 
    140 
    141        
    142     def __len__(self):
    143         return self.number_of_elements
    144 
    145     def __repr__(self):
    146         return 'Mesh: %d triangles, %d elements'\
    147                %(self.coordinates.shape[0], len(self))
    148 
    149 
    150 
    151     #DSG I don't know if I need this
    152     def build_vertexlist(self):
    153         pass
    154    
    155     # maybe add this latter..   
    156     def check_integrity(self):
    157         pass
    158    
    159 
    160     def get_normals(self):
    161         """Return all normal vectors.
    162         Return normal vectors for all triangles as an Nx6 array
    163         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
    164         """
    165         return self.normals
    166 
    167 
    168     def get_normal(self, i, j):
    169         """Return normal vector j of the i'th triangle.
    170         Return value is the numeric array slice [x, y]
    171         """
    172         return self.normals[i, 2*j:2*j+2]     
    173    
    174 
    175            
    176     def get_vertex_coordinates(self):
    177         """Return all vertex coordinates.
    178         Return all vertex coordinates for all triangles as an Nx6 array
    179         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)       
    180         """
    181         return self.vertex_coordinates
    182 
    183 
    184     def get_vertex_coordinate(self, i, j):
    185         """Return coordinates for vertex j of the i'th triangle.
    186         Return value is the numeric array slice [x, y]
    187         """
    188         return self.vertex_coordinates[i, 2*j:2*j+2]     
    189 
    190 
    191     def compute_vertex_coordinates(self):
    192         """Return vertex coordinates for all triangles as an Nx6 array
    193         (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
    194         """
    195 
    196         #FIXME: Perhaps they should be ordered as in obj files??
    197         #See quantity.get_vertex_values
    198        
    199         from Numeric import zeros, Float
    200 
    201         N = self.number_of_elements
    202         vertex_coordinates = zeros((N, 6), Float)
    203 
    204         for i in range(N):
    205             for j in range(3):
    206                 k = self.triangles[i,j]  #Index of vertex 0
    207                 v_k = self.coordinates[k]
    208                 vertex_coordinates[i, 2*j+0] = v_k[0]
    209                 vertex_coordinates[i, 2*j+1] = v_k[1]
    210 
    211         return vertex_coordinates
    212  
    213     def get_triangles(self, unique=False):
    214         """Get connectivity
    215         If unique is True give them only once as stored internally.
    216         If unique is False return
    217         """
    218 
    219         if unique is True:
    220             return self.triangles
    221         else:
    222             from Numeric import reshape, array, Int
    223                
    224             m = len(self)  #Number of volumes
    225             M = 3*m        #Total number of unique vertices
    226             return reshape(array(range(M)).astype(Int), (m,3))
    227 
    228 
    229 
     18from general_mesh import General_Mesh
    23019from Numeric import zeros, array, Float, Int, dot,transpose
    23120from LinearAlgebra import solve_linear_equations
    23221
    233 EPSILON = 1.0e-13
    23422
    235 class Interpolation(Mesh):
     23class Interpolation:
    23624
    237     def __init__(self, vertex_coordinates, triangles, point_coordinates):
     25    def __init__(self,
     26                 vertex_coordinates,
     27                 triangles,
     28                 point_coordinates = None):
     29
     30       
    23831        """ Build interpolation matrix mapping from
    23932        function values at vertices to function values at data points
    … …  
    25346
    25447        #Convert input to Numeric arrays
    255         point_coordinates = array(point_coordinates).astype(Float)
    25648        vertex_coordinates = array(vertex_coordinates).astype(Float)
    25749        triangles = array(triangles).astype(Int)               
    25850       
    25951        #Build underlying mesh
    260         Mesh.__init__(self, vertex_coordinates, triangles)
     52        self.mesh = General_Mesh(vertex_coordinates, triangles)
     53
     54        self.A_m = zeros((1,1), Float) #Not initialised
     55
     56        if point_coordinates:
     57            self.build_A_m(point_coordinates)
    26158       
     59    def build_A_m(self, point_coordinates):
    26260
     61        #Convert input to Numeric arrays
     62        point_coordinates = array(point_coordinates).astype(Float)
     63       
    26364        #Build n x m interpolation matrix       
    264         m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex)
     65        m = self.mesh.coordinates.shape[0] #Number of basis functions (1/vertex)
    26566        n = point_coordinates.shape[0]               #Number of data points         
    26667
    … …  
    27273
    27374            x = point_coordinates[i]
    274             for k in range(len(self)):
     75            for k in range(len(self.mesh)):
    27576                #For each triangle (brute force)
    27677                #FIXME: Real algorithm should only visit relevant triangles
    27778
    27879                #Get the three vertex_points
    279                 xi0 = self.get_vertex_coordinate(k, 0)
    280                 xi1 = self.get_vertex_coordinate(k, 1)
    281                 xi2 = self.get_vertex_coordinate(k, 2)                 
     80                xi0 = self.mesh.get_vertex_coordinate(k, 0)
     81                xi1 = self.mesh.get_vertex_coordinate(k, 1)
     82                xi2 = self.mesh.get_vertex_coordinate(k, 2)                 
    28283
    28384                #Get the three normals
    284                 n0 = self.get_normal(k, 0)
    285                 n1 = self.get_normal(k, 1)
    286                 n2 = self.get_normal(k, 2)               
     85                n0 = self.mesh.get_normal(k, 0)
     86                n1 = self.mesh.get_normal(k, 1)
     87                n2 = self.mesh.get_normal(k, 2)               
    28788
    28889                #Compute interpolation
    … …  
    299100
    300101                    #Assign values to A_m
    301                     j = self.triangles[k,0] #Global vertex id
     102                    j = self.mesh.triangles[k,0] #Global vertex id
    302103                    self.A_m[i, j] = sigma0
    303104
    304                     j = self.triangles[k,1] #Global vertex id
     105                    j = self.mesh.triangles[k,1] #Global vertex id
    305106                    self.A_m[i, j] = sigma1
    306107
    307                     j = self.triangles[k,2] #Global vertex id
     108                    j = self.mesh.triangles[k,2] #Global vertex id
    308109                    self.A_m[i, j] = sigma2
     110   
    309111
    310     def smooth_to_mesh(self, z):
    311         """ Linearly interpolate point data values to the vertices.
    312        
     112    def smooth_att_to_mesh(self, z):
     113        """ Linearly interpolate many points with an attribute to the vertices.
     114        Pre Condition:
     115          A_m has been initialised
     116         
    313117        Inputs:
    314118          z: Vector of data at the point_coordinates.  One value per point.
    … …  
    325129        return f
    326130
    327     def smooth_to_points(self, f):
    328         """ Linearly interpolate point data values to the vertices.
     131    def smooth_atts_to_mesh(self, z_m):
     132        """ Linearly interpolate many points with attributes to the vertices.
     133        Pre Condition:
     134          A_m has been initialised
     135         
     136        Inputs:
     137          z_m: A Matrix of data at the point_coordinates.
     138               One attribute per colum.
     139         
     140        """
     141        #Convert input to Numeric arrays
     142        z_m = array(z_m).astype(Float)
     143
     144       
     145        #Build n x m interpolation matrix       
     146        m = self.mesh.coordinates.shape[0] #Number of vertices
     147        n = z_m.shape[1]               #Number of data points         
     148
     149        f_m = zeros((n,m), Float)
     150       
     151        for i in range(z_m.shape[1]):
     152            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
     153       
     154       
     155    def smooth_att_to_points(self, f):
     156        """ Linearly interpolate vertices with an attribute to the points.
     157        Pre Condition:
     158          A_m has been initialised
    329159       
    330160        Inputs:
    … …  
    333163        """
    334164        return dot(self.A_m,f)
     165
     166    def smooth_atts_to_points(self, f_m):
     167        """ Linearly interpolate vertices with attributes to the points.
     168        Pre Condition:
     169          A_m has been initialised
     170       
     171        Inputs:
     172          f_m: Matrix of data at the  vertices.
     173               One attribute per colum.
     174         
     175        """
     176        #Convert input to Numeric arrays
     177        z_m = array(z_m).astype(Float)
     178
     179       
     180        #Build n x m interpolation matrix       
     181        m = self.mesh.coordinates.shape[0] #Number of vertices
     182        n = z_m.shape[1]               #Number of data points         
     183
     184        f_m = zeros((n,m), Float)
     185       
     186        for i in range(z_m.shape[1]):
     187            f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
     188       

  • inundation/ga/storm_surge/pyvolution/test_least_squares.py

    r316 r321  
    55
    66from least_squares import *
    7 from config import epsilon
     7from least_squares import EPSILON
    88from Numeric import allclose, array
    99
    … …  
    2828        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
    2929       
    30         mesh = Interpolation(points, vertices, data)
    31         assert allclose(mesh.A_m, [[1./3, 1./3, 1./3]])
     30        interp = Interpolation(points, vertices, data)
     31        assert allclose(interp.A_m, [[1./3, 1./3, 1./3]])
    3232       
    3333
    … …  
    4545        data = points #Use data at vertices       
    4646       
    47         mesh = Interpolation(points, vertices, data)
    48         assert allclose(mesh.A_m, [[1., 0., 0.],
     47        interp = Interpolation(points, vertices, data)
     48        assert allclose(interp.A_m, [[1., 0., 0.],
    4949                                      [0., 1., 0.],
    5050                                      [0., 0., 1.]])
    … …  
    6565        data = [ [0., 1.], [1., 0.], [1., 1.] ]
    6666       
    67         mesh = Interpolation(points, vertices, data)
    68        
    69         assert allclose(mesh.A_m, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
     67        interp = Interpolation(points, vertices, data)
     68       
     69        assert allclose(interp.A_m, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
    7070                                      [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
    7171                                      [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
    … …  
    8585        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
    8686       
    87         mesh = Interpolation(points, vertices, data)
    88        
    89         assert allclose(mesh.A_m, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0     
     87        interp = Interpolation(points, vertices, data)
     88       
     89        assert allclose(interp.A_m, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0     
    9090                                      [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
    9191                                      [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
    … …  
    105105        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
    106106       
    107         mesh = Interpolation(points, vertices, data)
    108 
    109         assert allclose(sum(mesh.A_m, axis=1), 1.0)
     107        interp = Interpolation(points, vertices, data)
     108
     109        assert allclose(sum(interp.A_m, axis=1), 1.0)
    110110       
    111111
    … …  
    125125        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
    126126
    127         mesh = Interpolation(points, triangles, data)
     127        interp = Interpolation(points, triangles, data)
    128128       
    129129        answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],  #Affects point d     
    … …  
    133133                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
    134134                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f       
    135         #print mesh.A_m
     135        #print interp.A_m
    136136        #print answer
    137137       
    138         assert allclose(mesh.A_m, answer)
    139 
    140     def test_smooth_to_mesh(self):
     138        assert allclose(interp.A_m, answer)
     139
     140    def test_smooth_att_to_mesh(self):
    141141        a = [0.0, 0.0]
    142142        b = [0.0, 5.0]
    … …  
    153153        data_coords = [ d1, d2, d3] #Use centroid as one data point       
    154154       
    155         mesh = Interpolation(points, triangles, data_coords)
     155        interp = Interpolation(points, triangles, data_coords)
    156156        z = [z1, z2, z3]
    157         f =  mesh.smooth_to_mesh(z)
     157        f =  interp.smooth_att_to_mesh(z)
    158158        answer = [0, 5., 5.]
    159159        assert allclose(f, answer)
    160160       
    161     def test_smooth_to_meshII(self):
     161    def test_smooth_att_to_meshII(self):
    162162        a = [0.0, 0.0]
    163163        b = [0.0, 5.0]
    … …  
    171171        data_coords = [ d1, d2, d3]       
    172172        z = self.linear_function(data_coords)
    173         mesh = Interpolation(points, triangles, data_coords)
    174         f =  mesh.smooth_to_mesh(z)
     173        interp = Interpolation(points, triangles, data_coords)
     174        f =  interp.smooth_att_to_mesh(z)
    175175        answer = self.linear_function(points)
    176176        assert allclose(f, answer)
    177177   
    178     def test_smooth_to_meshIII(self):
     178    def test_smooth_att_to_meshIII(self):
    179179        a = [-1.0, 0.0]
    180180        b = [3.0, 4.0]
    … …  
    202202       
    203203        z = self.linear_function(point_coords)
    204         mesh = Interpolation(vertices, triangles, point_coords)
    205         f =  mesh.smooth_to_mesh(z)
     204        interp = Interpolation(vertices, triangles, point_coords)
     205        f =  interp.smooth_att_to_mesh(z)
    206206        answer = self.linear_function(vertices)
    207207        assert allclose(f, answer)
    … …  
    211211        return point[:,0]+point[:,1]
    212212   
    213     def test_smooth_to_points(self):
     213    def test_smooth_att_to_points(self):
    214214        v0 = [0.0, 0.0]
    215215        v1 = [0.0, 5.0]
    … …  
    225225        point_coords = [ d0, d1, d2]
    226226       
    227         mesh = Interpolation(vertices, triangles, point_coords)
     227        interp = Interpolation(vertices, triangles, point_coords)
    228228        f = self.linear_function(vertices)
    229         z =  mesh.smooth_to_points(f)
     229        z =  interp.smooth_att_to_points(f)
    230230        answer = self.linear_function(point_coords)
    231231        #print "z",z
    … …  
    233233        assert allclose(z, answer)
    234234       
    235     def test_smooth_to_pointsII(self):
     235    def test_smooth_att_to_pointsII(self):
    236236        a = [-1.0, 0.0]
    237237        b = [3.0, 4.0]
    … …  
    258258                        [3.0,1.0]]
    259259       
    260         mesh = Interpolation(vertices, triangles, point_coords)
     260        interp = Interpolation(vertices, triangles, point_coords)
    261261        f = self.linear_function(vertices)
    262         z =  mesh.smooth_to_points(f)
     262        z =  interp.smooth_att_to_points(f)
    263263        answer = self.linear_function(point_coords)
    264264        #print "z",z
    265265        #print "answer",answer
    266         assert allclose(z, answer)
     266        assert allclose(z, answer)
     267
     268       
     269    def test_smooth_atts_to_mesh(self):
     270        a = [0.0, 0.0]
     271        b = [0.0, 5.0]
     272        c = [5.0, 0.0]
     273        points = [a, b, c]
     274        triangles = [ [1,0,2] ]   #bac
     275         
     276        d1 = [1.0, 1.0]
     277        d2 = [1.0, 3.0]
     278        d3 = [3.0,1.0]
     279        z1 = [2,4]
     280        z2 = [4,8]
     281        z3 = [4,8]
     282        data_coords = [ d1, d2, d3] #Use centroid as one data point       
     283       
     284        interp = Interpolation(points, triangles, data_coords)
     285        z = [z1, z2, z3]
     286        f =  interp.smooth_att_to_mesh(z)
     287        answer = [[0,0], [5., 10.], [5., 10.]]
     288        assert allclose(f, answer)
     289       
    267290#-------------------------------------------------------------
    268291if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.