Changeset 328


Ignore:
Timestamp:
Sep 21, 2004, 1:26:23 PM (20 years ago)
Author:
duncan
Message:

adding reading .tsh files to smoothing

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

Legend:

Unmodified
Added
Removed

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

    r321 r328  
    1515EPSILON = 1.0e-13
    1616
    17 
    1817from general_mesh import General_Mesh
    1918from Numeric import zeros, array, Float, Int, dot,transpose
    2019from LinearAlgebra import solve_linear_equations
     20
     21
     22def smooth_attributes_to_mesh(vertex_coordinates,
     23                              triangles,
     24                              point_coordinates,
     25                              point_attributes):
     26    """
     27    """
     28    interp = Interpolation(vertex_coordinates,
     29                              triangles,
     30                              point_coordinates)
     31    vertex_attributes = interp.smooth_attributes_to_mesh(point_attributes)
     32    return vertex_attributes
    2133
    2234
    … …  
    5264        self.mesh = General_Mesh(vertex_coordinates, triangles)
    5365
    54         self.A_m = zeros((1,1), Float) #Not initialised
     66        self.A = zeros((1,1), Float) #Not initialised
    5567
    5668        if point_coordinates:
    57             self.build_A_m(point_coordinates)
     69            self.build_A(point_coordinates)
    5870       
    59     def build_A_m(self, point_coordinates):
     71    def build_A(self, point_coordinates):
     72        """
     73        "A" is the matrix of hat functions.  It has 1 row per data point and 1 column per vertex.
     74
     75        Post Condition
     76        "A" has been initialised
     77        """
    6078
    6179        #Convert input to Numeric arrays
    … …  
    6684        n = point_coordinates.shape[0]               #Number of data points         
    6785
    68         self.A_m = zeros((n,m), Float)
     86        self.A = zeros((n,m), Float)
    6987
    7088        #Compute matrix elements
    … …  
    99117                if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
    100118
    101                     #Assign values to A_m
     119                    #Assign values to A
    102120                    j = self.mesh.triangles[k,0] #Global vertex id
    103                     self.A_m[i, j] = sigma0
     121                    self.A[i, j] = sigma0
    104122
    105123                    j = self.mesh.triangles[k,1] #Global vertex id
    106                     self.A_m[i, j] = sigma1
     124                    self.A[i, j] = sigma1
    107125
    108126                    j = self.mesh.triangles[k,2] #Global vertex id
    109                     self.A_m[i, j] = sigma2
     127                    self.A[i, j] = sigma2
    110128   
    111129
    112     def smooth_att_to_mesh(self, z):
     130    def smooth_attributes_to_mesh(self, z):
    113131        """ Linearly interpolate many points with an attribute to the vertices.
    114132        Pre Condition:
    115           A_m has been initialised
     133          "A" has been initialised
    116134         
    117135        Inputs:
    … …  
    121139        #Convert input to Numeric arrays
    122140        z = array(z).astype(Float)
    123         At_m = transpose(self.A_m)
     141        At = transpose(self.A)
    124142        #print "z", z
    125         #print "At_m",At_m
    126         self.AtA_m = dot(At_m, self.A_m)
    127         Atz_m = dot(At_m, z)
    128         f = solve_linear_equations(self.AtA_m,Atz_m)
     143        #print "At",At
     144        self.AtA = dot(At, self.A)
     145        Atz = dot(At, z)
     146        f = solve_linear_equations(self.AtA,Atz)
    129147        return f
    130148
    131     def smooth_atts_to_mesh(self, z_m):
    132         """ Linearly interpolate many points with attributes to the vertices.
     149       
     150    def interpolate_attributes_to_points(self, f):
     151        """ Linearly interpolate vertices with attributes to the points.
    133152        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
     153          A has been initialised
    159154       
    160155        Inputs:
    161           z: Vector of data at the point_coordinates.  One value per point.
    162          
    163         """
    164         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.
     156          f: Matrix of data at the  vertices.
    173157               One attribute per colum.
    174158         
    175159        """
    176         #Convert input to Numeric arrays
    177         z_m = array(z_m).astype(Float)
     160        return dot(self.A,f)
    178161
     162#-------------------------------------------------------------
     163if __name__ == "__main__":
     164   
     165    import os, sys
     166    from load_mesh.loadASCII import mesh_file_to_mesh_dictionary
     167    usage = "usage: %s pmesh_file_name" %         os.path.basename(sys.argv[0])
     168
     169    if len(sys.argv) < 2:
     170        print usage
     171    else:
     172        mesh_file_name = sys.argv[1]
     173
     174        meshdic = mesh_file_to_mesh_dictionary(mesh_file_name)
     175        vertex_coordinates = meshdic['generatedpointlist']
     176        triangles = meshdic['generatedtrianglelist']
    179177       
    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)
     178        d1 = [1.0, 1.0]
     179        d2 = [1.0, 3.0]
     180        d3 = [3.0,1.0]
     181        z1 = 2
     182        z2 = 4
     183        z3 = 4
     184        data_coords = [ d1, d2, d3]
     185        z = [z1, z2, z3]
    185186       
    186         for i in range(z_m.shape[1]):
    187             f_m[:,i] = smooth_att_to_mesh(z_m[:,i])
     187        f = smooth_attributes_to_mesh(vertex_coordinates,
     188                                      triangles,
     189                                      data_coords,
     190                                      z)
     191        print f
    188192       
     193       

  • inundation/ga/storm_surge/pyvolution/load_mesh/loadASCII.py

    r298 r328  
    11
    22from string import  find, rfind
    3    
     3
     4def mesh_file_to_mesh_dictionary(fileName):
     5    """Load a pmesh file.  Returning the mesh dictionary.
     6    """
     7    try:
     8        meshdic = import_trianglulation(fileName)
     9    except IOError, e:       
     10        msg = 'Could not open file %s ' %fileName
     11        raise IOError, msg
     12    return meshdic
     13
    414def import_trianglulation(ofile):
    515    """

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

    r298 r328  
    11def pmesh_to_domain(fileName, tags = None, setting_function = None):
     2    """
     3    """
     4    #FIXME:  The current design doesn't seem to accomadate tags and
     5    # setting_functions being passed into the domain at this point.
     6
     7    #FIXME: Plus the names of the functions are no longer appropriate,
     8    #since domain objects aren't being returned.
     9   
    210    import sys
    311    from config import pmesh_filename
    … …  
    1220
    1321       
    14     domain  = pmesh_dictionary_to_domain(meshdic,
    15                                          setting_function = setting_function)
    16     if tags:
    17         domain.set_boundary(tags)
    18        
    19     return domain
     22    return pmesh_dictionary_to_domain(meshdic,
     23                                      setting_function = setting_function)
    2024
    2125

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

    r321 r328  
    55
    66from least_squares import *
    7 from least_squares import EPSILON
    87from Numeric import allclose, array
    98
    … …  
    2928       
    3029        interp = Interpolation(points, vertices, data)
    31         assert allclose(interp.A_m, [[1./3, 1./3, 1./3]])
     30        assert allclose(interp.A, [[1./3, 1./3, 1./3]])
    3231       
    3332
    … …  
    4645       
    4746        interp = Interpolation(points, vertices, data)
    48         assert allclose(interp.A_m, [[1., 0., 0.],
     47        assert allclose(interp.A, [[1., 0., 0.],
    4948                                      [0., 1., 0.],
    5049                                      [0., 0., 1.]])
    … …  
    6766        interp = Interpolation(points, vertices, data)
    6867       
    69         assert allclose(interp.A_m, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
     68        assert allclose(interp.A, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
    7069                                      [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
    7170                                      [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2
    … …  
    8786        interp = Interpolation(points, vertices, data)
    8887       
    89         assert allclose(interp.A_m, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0     
     88        assert allclose(interp.A, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0     
    9089                                      [0.25, 0.0, 0.75],  #Affects vertex 0 and 2
    9190                                      [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2
    … …  
    107106        interp = Interpolation(points, vertices, data)
    108107
    109         assert allclose(sum(interp.A_m, axis=1), 1.0)
     108        assert allclose(sum(interp.A, axis=1), 1.0)
    110109       
    111110
    … …  
    133132                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
    134133                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f       
    135         #print interp.A_m
     134        #print interp.A
    136135        #print answer
    137136       
    138         assert allclose(interp.A_m, answer)
    139 
    140     def test_smooth_att_to_mesh(self):
     137        assert allclose(interp.A, answer)
     138
     139    def test_smooth_attributes_to_mesh(self):
    141140        a = [0.0, 0.0]
    142141        b = [0.0, 5.0]
    … …  
    151150        z2 = 4
    152151        z3 = 4
    153         data_coords = [ d1, d2, d3] #Use centroid as one data point       
     152        data_coords = [ d1, d2, d3]
    154153       
    155154        interp = Interpolation(points, triangles, data_coords)
    156155        z = [z1, z2, z3]
    157         f =  interp.smooth_att_to_mesh(z)
     156        f =  interp.smooth_attributes_to_mesh(z)
    158157        answer = [0, 5., 5.]
    159158        assert allclose(f, answer)
    160159       
    161     def test_smooth_att_to_meshII(self):
     160    def test_smooth_attributes_to_meshII(self):
    162161        a = [0.0, 0.0]
    163162        b = [0.0, 5.0]
    … …  
    172171        z = self.linear_function(data_coords)
    173172        interp = Interpolation(points, triangles, data_coords)
    174         f =  interp.smooth_att_to_mesh(z)
     173        f =  interp.smooth_attributes_to_mesh(z)
    175174        answer = self.linear_function(points)
    176175        assert allclose(f, answer)
    177176   
    178     def test_smooth_att_to_meshIII(self):
     177    def test_smooth_attributes_to_meshIII(self):
    179178        a = [-1.0, 0.0]
    180179        b = [3.0, 4.0]
    … …  
    203202        z = self.linear_function(point_coords)
    204203        interp = Interpolation(vertices, triangles, point_coords)
    205         f =  interp.smooth_att_to_mesh(z)
     204        f =  interp.smooth_attributes_to_mesh(z)
    206205        answer = self.linear_function(vertices)
    207206        assert allclose(f, answer)
    … …  
    210209        point = array(point)
    211210        return point[:,0]+point[:,1]
    212    
    213     def test_smooth_att_to_points(self):
     211
     212    def test_smooth_attributes_to_meshIV(self):
     213        """ Testing 2 attributes smoothed to the mesh
     214        """
     215        a = [0.0, 0.0]
     216        b = [0.0, 5.0]
     217        c = [5.0, 0.0]
     218        points = [a, b, c]
     219        triangles = [ [1,0,2] ]   #bac
     220         
     221        d1 = [1.0, 1.0]
     222        d2 = [1.0, 3.0]
     223        d3 = [3.0,1.0]
     224        z1 = [2,4]
     225        z2 = [4,8]
     226        z3 = [4,8]
     227        data_coords = [ d1, d2, d3]
     228       
     229        interp = Interpolation(points, triangles, data_coords)
     230        z = [z1, z2, z3]
     231        f =  interp.smooth_attributes_to_mesh(z)
     232        answer = [[0,0], [5., 10.], [5., 10.]]
     233        assert allclose(f, answer)
     234       
     235    def test_interpolate_attributes_to_points(self):
    214236        v0 = [0.0, 0.0]
    215237        v1 = [0.0, 5.0]
    … …  
    227249        interp = Interpolation(vertices, triangles, point_coords)
    228250        f = self.linear_function(vertices)
    229         z =  interp.smooth_att_to_points(f)
     251        z =  interp.interpolate_attributes_to_points(f)
    230252        answer = self.linear_function(point_coords)
    231253        #print "z",z
    … …  
    233255        assert allclose(z, answer)
    234256       
    235     def test_smooth_att_to_pointsII(self):
     257    def test_interpolate_attributes_to_pointsII(self):
    236258        a = [-1.0, 0.0]
    237259        b = [3.0, 4.0]
    … …  
    260282        interp = Interpolation(vertices, triangles, point_coords)
    261283        f = self.linear_function(vertices)
    262         z =  interp.smooth_att_to_points(f)
     284        z =  interp.interpolate_attributes_to_points(f)
    263285        answer = self.linear_function(point_coords)
    264286        #print "z",z
    … …  
    267289
    268290       
    269     def test_smooth_atts_to_mesh(self):
     291       
     292    def test_smooth_attributes_to_mesh_function(self):
     293        """ Testing 2 attributes smoothed to the mesh
     294        """
    270295        a = [0.0, 0.0]
    271296        b = [0.0, 5.0]
    … …  
    280305        z2 = [4,8]
    281306        z3 = [4,8]
    282         data_coords = [ d1, d2, d3] #Use centroid as one data point       
    283        
    284         interp = Interpolation(points, triangles, data_coords)
     307        data_coords = [ d1, d2, d3]
     308       
    285309        z = [z1, z2, z3]
    286         f =  interp.smooth_att_to_mesh(z)
     310        f =  smooth_attributes_to_mesh(points, triangles, data_coords,z)
    287311        answer = [[0,0], [5., 10.], [5., 10.]]
    288312        assert allclose(f, answer)
Note: See TracChangeset for help on using the changeset viewer.