Changeset 328
- Timestamp:
- Sep 21, 2004, 1:26:23 PM (20 years ago)
- 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 15 15 EPSILON = 1.0e-13 16 16 17 18 17 from general_mesh import General_Mesh 19 18 from Numeric import zeros, array, Float, Int, dot,transpose 20 19 from LinearAlgebra import solve_linear_equations 20 21 22 def 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 21 33 22 34 … … 52 64 self.mesh = General_Mesh(vertex_coordinates, triangles) 53 65 54 self.A _m= zeros((1,1), Float) #Not initialised66 self.A = zeros((1,1), Float) #Not initialised 55 67 56 68 if point_coordinates: 57 self.build_A _m(point_coordinates)69 self.build_A(point_coordinates) 58 70 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 """ 60 78 61 79 #Convert input to Numeric arrays … … 66 84 n = point_coordinates.shape[0] #Number of data points 67 85 68 self.A _m= zeros((n,m), Float)86 self.A = zeros((n,m), Float) 69 87 70 88 #Compute matrix elements … … 99 117 if sigma0 >= 0 and sigma1 >= 0 and sigma2 >= 0: 100 118 101 #Assign values to A _m119 #Assign values to A 102 120 j = self.mesh.triangles[k,0] #Global vertex id 103 self.A _m[i, j] = sigma0121 self.A[i, j] = sigma0 104 122 105 123 j = self.mesh.triangles[k,1] #Global vertex id 106 self.A _m[i, j] = sigma1124 self.A[i, j] = sigma1 107 125 108 126 j = self.mesh.triangles[k,2] #Global vertex id 109 self.A _m[i, j] = sigma2127 self.A[i, j] = sigma2 110 128 111 129 112 def smooth_att _to_mesh(self, z):130 def smooth_attributes_to_mesh(self, z): 113 131 """ Linearly interpolate many points with an attribute to the vertices. 114 132 Pre Condition: 115 A_mhas been initialised133 "A" has been initialised 116 134 117 135 Inputs: … … 121 139 #Convert input to Numeric arrays 122 140 z = array(z).astype(Float) 123 At _m = transpose(self.A_m)141 At = transpose(self.A) 124 142 #print "z", z 125 #print "At _m",At_m126 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) 129 147 return f 130 148 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. 133 152 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 159 154 160 155 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. 173 157 One attribute per colum. 174 158 175 159 """ 176 #Convert input to Numeric arrays 177 z_m = array(z_m).astype(Float) 160 return dot(self.A,f) 178 161 162 #------------------------------------------------------------- 163 if __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'] 179 177 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] 185 186 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 188 192 193 -
inundation/ga/storm_surge/pyvolution/load_mesh/loadASCII.py
r298 r328 1 1 2 2 from string import find, rfind 3 3 4 def 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 4 14 def import_trianglulation(ofile): 5 15 """ -
inundation/ga/storm_surge/pyvolution/pmesh2domain.py
r298 r328 1 1 def 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 2 10 import sys 3 11 from config import pmesh_filename … … 12 20 13 21 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) 20 24 21 25 -
inundation/ga/storm_surge/pyvolution/test_least_squares.py
r321 r328 5 5 6 6 from least_squares import * 7 from least_squares import EPSILON8 7 from Numeric import allclose, array 9 8 … … 29 28 30 29 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]]) 32 31 33 32 … … 46 45 47 46 interp = Interpolation(points, vertices, data) 48 assert allclose(interp.A _m, [[1., 0., 0.],47 assert allclose(interp.A, [[1., 0., 0.], 49 48 [0., 1., 0.], 50 49 [0., 0., 1.]]) … … 67 66 interp = Interpolation(points, vertices, data) 68 67 69 assert allclose(interp.A _m, [[0.5, 0.5, 0.0], #Affects vertex 1 and 068 assert allclose(interp.A, [[0.5, 0.5, 0.0], #Affects vertex 1 and 0 70 69 [0.5, 0.0, 0.5], #Affects vertex 0 and 2 71 70 [0.0, 0.5, 0.5]]) #Affects vertex 1 and 2 … … 87 86 interp = Interpolation(points, vertices, data) 88 87 89 assert allclose(interp.A _m, [[0.25, 0.75, 0.0], #Affects vertex 1 and 088 assert allclose(interp.A, [[0.25, 0.75, 0.0], #Affects vertex 1 and 0 90 89 [0.25, 0.0, 0.75], #Affects vertex 0 and 2 91 90 [0.0, 0.25, 0.75]]) #Affects vertex 1 and 2 … … 107 106 interp = Interpolation(points, vertices, data) 108 107 109 assert allclose(sum(interp.A _m, axis=1), 1.0)108 assert allclose(sum(interp.A, axis=1), 1.0) 110 109 111 110 … … 133 132 [0.25, 0.75, 0.0, 0.0, 0.0, 0.0], #Affects points a and b 134 133 [1./3, 0.0, 0.0, 0.0, 1./3, 1./3]] #Affects points a,e and f 135 #print interp.A _m134 #print interp.A 136 135 #print answer 137 136 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): 141 140 a = [0.0, 0.0] 142 141 b = [0.0, 5.0] … … 151 150 z2 = 4 152 151 z3 = 4 153 data_coords = [ d1, d2, d3] #Use centroid as one data point152 data_coords = [ d1, d2, d3] 154 153 155 154 interp = Interpolation(points, triangles, data_coords) 156 155 z = [z1, z2, z3] 157 f = interp.smooth_att _to_mesh(z)156 f = interp.smooth_attributes_to_mesh(z) 158 157 answer = [0, 5., 5.] 159 158 assert allclose(f, answer) 160 159 161 def test_smooth_att _to_meshII(self):160 def test_smooth_attributes_to_meshII(self): 162 161 a = [0.0, 0.0] 163 162 b = [0.0, 5.0] … … 172 171 z = self.linear_function(data_coords) 173 172 interp = Interpolation(points, triangles, data_coords) 174 f = interp.smooth_att _to_mesh(z)173 f = interp.smooth_attributes_to_mesh(z) 175 174 answer = self.linear_function(points) 176 175 assert allclose(f, answer) 177 176 178 def test_smooth_att _to_meshIII(self):177 def test_smooth_attributes_to_meshIII(self): 179 178 a = [-1.0, 0.0] 180 179 b = [3.0, 4.0] … … 203 202 z = self.linear_function(point_coords) 204 203 interp = Interpolation(vertices, triangles, point_coords) 205 f = interp.smooth_att _to_mesh(z)204 f = interp.smooth_attributes_to_mesh(z) 206 205 answer = self.linear_function(vertices) 207 206 assert allclose(f, answer) … … 210 209 point = array(point) 211 210 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): 214 236 v0 = [0.0, 0.0] 215 237 v1 = [0.0, 5.0] … … 227 249 interp = Interpolation(vertices, triangles, point_coords) 228 250 f = self.linear_function(vertices) 229 z = interp. smooth_att_to_points(f)251 z = interp.interpolate_attributes_to_points(f) 230 252 answer = self.linear_function(point_coords) 231 253 #print "z",z … … 233 255 assert allclose(z, answer) 234 256 235 def test_ smooth_att_to_pointsII(self):257 def test_interpolate_attributes_to_pointsII(self): 236 258 a = [-1.0, 0.0] 237 259 b = [3.0, 4.0] … … 260 282 interp = Interpolation(vertices, triangles, point_coords) 261 283 f = self.linear_function(vertices) 262 z = interp. smooth_att_to_points(f)284 z = interp.interpolate_attributes_to_points(f) 263 285 answer = self.linear_function(point_coords) 264 286 #print "z",z … … 267 289 268 290 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 """ 270 295 a = [0.0, 0.0] 271 296 b = [0.0, 5.0] … … 280 305 z2 = [4,8] 281 306 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 285 309 z = [z1, z2, z3] 286 f = interp.smooth_att_to_mesh(z)310 f = smooth_attributes_to_mesh(points, triangles, data_coords,z) 287 311 answer = [[0,0], [5., 10.], [5., 10.]] 288 312 assert allclose(f, answer)
Note: See TracChangeset
for help on using the changeset viewer.