Changeset 309
- Timestamp:
- Sep 16, 2004, 12:01:21 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/least_squares.py
r305 r309 12 12 class Interpolation(Mesh): 13 13 14 def __init__(self, vertex_coordinates, triangles, data ):14 def __init__(self, vertex_coordinates, triangles, data_coordinates): 15 15 """ Build interpolation matrix mapping from 16 16 function values at vertices to function values at data points … … 24 24 integers representing indices of all vertices in the mesh. 25 25 26 data : List of coordinate pairs [x, y] of data points26 data_coordinates: List of coordinate pairs [x, y] of data points 27 27 (or an nx2 Numeric array) 28 28 … … 33 33 34 34 #Convert input to Numeric arrays 35 data = array(data).astype(Float)35 data_coordinates = array(data_coordinates).astype(Float) 36 36 vertex_coordinates = array(vertex_coordinates).astype(Float) 37 37 triangles = array(triangles).astype(Int) … … 43 43 #Build n x m interpolation matrix 44 44 m = vertex_coordinates.shape[0] #Number of basis functions (1/vertex) 45 n = data .shape[0] #Number of data points45 n = data_coordinates.shape[0] #Number of data points 46 46 47 self. matrix= zeros((n,m), Float)47 self.A_m = zeros((n,m), Float) 48 48 49 49 #Compute matrix elements 50 50 for i in range(n): 51 #For each data point51 #For each data_coordinate point 52 52 53 x = data [i]53 x = data_coordinates[i] 54 54 for k in range(len(self)): 55 55 #For each triangle (brute force) … … 78 78 if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0: 79 79 80 #Assign values to matrix80 #Assign values to A_m 81 81 j = self.triangles[k,0] #Global vertex id 82 self. matrix[i, j] = sigma082 self.A_m[i, j] = sigma0 83 83 84 84 j = self.triangles[k,1] #Global vertex id 85 self. matrix[i, j] = sigma185 self.A_m[i, j] = sigma1 86 86 87 87 j = self.triangles[k,2] #Global vertex id 88 self. matrix[i, j] = sigma288 self.A_m[i, j] = sigma2 89 89 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 38 38 domain.filename = 'weir' 39 39 domain.checkpoint = False #True 40 domain.visualise = False#True #False40 domain.visualise = True #False 41 41 domain.default_order = 2 42 42 -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r305 r309 24 24 c = [2.0,0.0] 25 25 points = [a, b, c] 26 triangles = [ [1,0,2] ] #bac26 vertices = [ [1,0,2] ] #bac 27 27 28 28 data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point 29 29 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]]) 32 32 33 33 … … 41 41 c = [2.0,0.0] 42 42 points = [a, b, c] 43 triangles = [ [1,0,2] ] #bac43 vertices = [ [1,0,2] ] #bac 44 44 45 45 data = points #Use data at vertices 46 46 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.], 49 49 [0., 1., 0.], 50 50 [0., 0., 1.]]) … … 61 61 c = [2.0,0.0] 62 62 points = [a, b, c] 63 triangles = [ [1,0,2] ] #bac63 vertices = [ [1,0,2] ] #bac 64 64 65 65 data = [ [0., 1.], [1., 0.], [1., 1.] ] 66 66 67 mesh = Interpolation(points, triangles, data)67 mesh = Interpolation(points, vertices, data) 68 68 69 assert allclose(mesh. matrix, [[0.5, 0.5, 0.0], #Affects vertex 1 and 069 assert allclose(mesh.A_m, [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 70 70 [0.5, 0.0, 0.5], #Affects vertex 0 and 2 71 71 [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2 … … 81 81 c = [2.0,0.0] 82 82 points = [a, b, c] 83 triangles = [ [1,0,2] ] #bac83 vertices = [ [1,0,2] ] #bac 84 84 85 85 data = [ [0., 1.5], [1.5, 0.], [1.5, 0.5] ] 86 86 87 mesh = Interpolation(points, triangles, data)87 mesh = Interpolation(points, vertices, data) 88 88 89 assert allclose(mesh. matrix, [[0.25, 0.75, 0.0], #Affects vertex 1 and 089 assert allclose(mesh.A_m, [[0.25, 0.75, 0.0], #Affects vertex 1 and 0 90 90 [0.25, 0.0, 0.75], #Affects vertex 0 and 2 91 91 [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2 … … 101 101 c = [2.0,0.0] 102 102 points = [a, b, c] 103 triangles = [ [1,0,2] ] #bac103 vertices = [ [1,0,2] ] #bac 104 104 105 105 data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ] 106 106 107 mesh = Interpolation(points, triangles, data)107 mesh = Interpolation(points, vertices, data) 108 108 109 assert allclose(sum(mesh. matrix, axis=1), 1.0)109 assert allclose(sum(mesh.A_m, axis=1), 1.0) 110 110 111 111 112 112 113 113 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 114 120 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 128 123 129 124 #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] ] 131 126 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) 133 139 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 137 154 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) 138 160 139 161 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.