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.