Changeset 309


Ignore:
Timestamp:
Sep 16, 2004, 12:01:21 PM (21 years ago)
Author:
duncan
Message:

added basic smoothing function to least squares

Location:
inundation/ga/storm_surge/pyvolution
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r305 r309  
    1212class Interpolation(Mesh):
    1313
    14     def __init__(self, vertex_coordinates, triangles, data):
     14    def __init__(self, vertex_coordinates, triangles, data_coordinates):
    1515        """ Build interpolation matrix mapping from
    1616        function values at vertices to function values at data points
    … …  
    2424          integers representing indices of all vertices in the mesh.
    2525
    26           data: List of coordinate pairs [x, y] of data points
     26          data_coordinates: List of coordinate pairs [x, y] of data points
    2727          (or an nx2 Numeric array)
    2828         
    … …  
    3333
    3434        #Convert input to Numeric arrays
    35         data = array(data).astype(Float)
     35        data_coordinates = array(data_coordinates).astype(Float)
    3636        vertex_coordinates = array(vertex_coordinates).astype(Float)
    3737        triangles = array(triangles).astype(Int)               
    … …  
    4343        #Build n x m interpolation matrix       
    4444        m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex)
    45         n = data.shape[0]               #Number of data points         
     45        n = data_coordinates.shape[0]               #Number of data points         
    4646
    47         self.matrix = zeros((n,m), Float)
     47        self.A_m = zeros((n,m), Float)
    4848
    4949        #Compute matrix elements
    5050        for i in range(n):
    51             #For each data point
     51            #For each data_coordinate point
    5252
    53             x = data[i]
     53            x = data_coordinates[i]
    5454            for k in range(len(self)):
    5555                #For each triangle (brute force)
    … …  
    7878                if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0:
    7979
    80                     #Assign values to matrix
     80                    #Assign values to A_m
    8181                    j = self.triangles[k,0] #Global vertex id
    82                     self.matrix[i, j] = sigma0
     82                    self.A_m[i, j] = sigma0
    8383
    8484                    j = self.triangles[k,1] #Global vertex id
    85                     self.matrix[i, j] = sigma1
     85                    self.A_m[i, j] = sigma1
    8686
    8787                    j = self.triangles[k,2] #Global vertex id
    88                     self.matrix[i, j] = sigma2
     88                    self.A_m[i, j] = sigma2
    8989
    90                
    91 
    92                
    93 
    94                
    95                
    96                
     90    def smooth_to_mesh(self, z):
     91        from Numeric import zeros, array, Float,transpose,dot
     92        from LinearAlgebra import solve_linear_equations
     93        #Convert input to Numeric arrays
     94        z = array(z).astype(Float)
     95        At_m = transpose(self.A_m)
     96        #print "z", z
     97        #print "At_m",At_m
     98        self.AtA_m = dot(At_m, self.A_m)
     99        Atz_m = dot(At_m, z)
     100        f = solve_linear_equations(self.AtA_m,Atz_m)
     101        return f

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

    r300 r309  
    3838domain.filename = 'weir'
    3939domain.checkpoint = False #True
    40 domain.visualise = False#True #False
     40domain.visualise = True #False
    4141domain.default_order = 2
    4242

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

    r305 r309  
    2424        c = [2.0,0.0]
    2525        points = [a, b, c]
    26         triangles = [ [1,0,2] ]   #bac
     26        vertices = [ [1,0,2] ]   #bac
    2727
    2828        data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point       
    2929       
    30         mesh = Interpolation(points, triangles, data)
    31         assert allclose(mesh.matrix, [[1./3, 1./3, 1./3]])
     30        mesh = Interpolation(points, vertices, data)
     31        assert allclose(mesh.A_m, [[1./3, 1./3, 1./3]])
    3232       
    3333
    … …  
    4141        c = [2.0,0.0]
    4242        points = [a, b, c]
    43         triangles = [ [1,0,2] ]   #bac
     43        vertices = [ [1,0,2] ]   #bac
    4444
    4545        data = points #Use data at vertices       
    4646       
    47         mesh = Interpolation(points, triangles, data)
    48         assert allclose(mesh.matrix, [[1., 0., 0.],
     47        mesh = Interpolation(points, vertices, data)
     48        assert allclose(mesh.A_m, [[1., 0., 0.],
    4949                                      [0., 1., 0.],
    5050                                      [0., 0., 1.]])
    … …  
    6161        c = [2.0,0.0]
    6262        points = [a, b, c]
    63         triangles = [ [1,0,2] ]   #bac
     63        vertices = [ [1,0,2] ]   #bac
    6464
    6565        data = [ [0., 1.], [1., 0.], [1., 1.] ]
    6666       
    67         mesh = Interpolation(points, triangles, data)
     67        mesh = Interpolation(points, vertices, data)
    6868       
    69         assert allclose(mesh.matrix, [[0.5, 0.5, 0.0],  #Affects vertex 1 and 0       
     69        assert allclose(mesh.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
    … …  
    8181        c = [2.0,0.0]
    8282        points = [a, b, c]
    83         triangles = [ [1,0,2] ]   #bac
     83        vertices = [ [1,0,2] ]   #bac
    8484
    8585        data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ]
    8686       
    87         mesh = Interpolation(points, triangles, data)
     87        mesh = Interpolation(points, vertices, data)
    8888       
    89         assert allclose(mesh.matrix, [[0.25, 0.75, 0.0],  #Affects vertex 1 and 0     
     89        assert allclose(mesh.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
    … …  
    101101        c = [2.0,0.0]
    102102        points = [a, b, c]
    103         triangles = [ [1,0,2] ]   #bac
     103        vertices = [ [1,0,2] ]   #bac
    104104
    105105        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
    106106       
    107         mesh = Interpolation(points, triangles, data)
     107        mesh = Interpolation(points, vertices, data)
    108108
    109         assert allclose(sum(mesh.matrix, axis=1), 1.0)
     109        assert allclose(sum(mesh.A_m, axis=1), 1.0)
    110110       
    111111
    112112           
    113113    def test_more_triangles(self):
     114        a = [-1.0, 0.0]
     115        b = [3.0, 4.0]
     116        c = [4.0,1.0]
     117        d = [-3.0, 2.0] #3
     118        e = [-1.0,-2.0]
     119        f = [1.0, -2.0] #5
    114120
    115 
    116         #FIXME: Not finished yet
    117         a = [0.0, 0.0]
    118         b = [0.0, 2.0]
    119         c = [2.0,0.0]
    120         d = [0.0, 4.0]
    121         e = [2.0, 2.0]
    122         f = [4.0,0.0]
    123 
    124         points = [a, b, c, d, e, f]
    125         #bac, bce, ecf, dbe, daf, dae
    126         triangles = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]]
    127 
     121        points = [a, b, c, d,e,f]
     122        triangles = [[0,1,3],[1,0,2],[0,4,5], [0,5,2]] #abd bac aef afc
    128123
    129124        #Data points
    130         data = [ [0., 2.0], [1.15, 0.75], [0.8888, 0.21], [1.000001, 0.00001] ]
     125        data = [ [-3., 2.0], [-2, 1], [0.0, 1], [0, 3], [2, 3], [-1.0/3,-4./3] ]
    131126
    132         mesh = Interpolation(points, triangles, data)       
     127        mesh = Interpolation(points, triangles, data)
     128       
     129        answer = [[0.0, 0.0, 0.0,  1.0, 0.0, 0.0],  #Affects point d     
     130                  [0.5, 0.0, 0.0, 0.5, 0.0, 0.0],  #Affects points a and d
     131                  [0.75, 0.25, 0.0, 0.0, 0.0, 0.0],  #Affects points a and b
     132                  [0.0, 0.5, 0.0, 0.5, 0.0, 0.0],  #Affects points a and d
     133                  [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b
     134                  [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f       
     135        #print mesh.A_m
     136        #print answer
     137       
     138        assert allclose(mesh.A_m, answer)
    133139
    134         #print mesh.matrix
    135        
    136 
     140    def test_smooth_to_mesh(self):
     141        a = [0.0, 0.0]
     142        b = [0.0, 5.0]
     143        c = [5.0, 0.0]
     144        points = [a, b, c]
     145        triangles = [ [1,0,2] ]   #bac
     146         
     147        d1 = [1.0, 1.0]
     148        d2 = [1.0, 3.0]
     149        d3 = [3.0,1.0]
     150        z1 = 2
     151        z2 = 4
     152        z3 = 4
     153        data_coords = [ d1, d2, d3] #Use centroid as one data point       
    137154       
     155        mesh = Interpolation(points, triangles, data_coords)
     156        z = [z1, z2, z3]
     157        f =  mesh.smooth_to_mesh(z)
     158        answer = [0, 5., 5.]
     159        assert allclose(f, answer)
    138160       
    139161#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.